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1.  INTRODUCTION 


High  velocity  impact  loading  generates  high  pressures,  very  high  strain 
rates,  and  high  temperatures.  Most  structural  alloys  exhibit  strain  rate  de¬ 
pendent  material  behavior  under  high  strain  rates  (>1000/sec) .  The 
constitutive  description  of  structural  alloys  under  impact  loading  must  there¬ 
fore  incorporate  the  effects  of  pressure,  strain  rate,  temperature,  and 
loading  history  on  strength  and  plastic  deformation.  Material  strength  in¬ 
creases  with  increasing  strain  rates.  Depending  on  the  various  levels  and 
conditions  of  these  variables,  the  failure  mechanisms  can  change  widely.  For 
instance,  at  rapid  plastic  deformations,  the  thermal  softening  of  materials 
enhances  the  formation  of  localized  adiabatic  shear  bands.  High  tensile  pres¬ 
sures  are  often  generated  due  to  wave  interactions.  These  pressures  can  cause 
void  nucleation  and  growth.  The  use  of  advanced,  continuum  mechanics  based, 
three  dimensional,  thermoviscoplastic  constitutive  and  failure  models  incor¬ 
porating  these  physical  features  is  essential  for  accurate  computer  code 
calculations  of  practical  impact  engineering  problems. 

Ductile  material  failure  is  often  accompanied  by  nucleation  and  growth 
of  microvoids.  Most  initially  intact  structural  alloys  develop 
microvoids/microcracks  during  plastic  deformation  under  both  quasi-static  and 
dynamic  loading  conditions.  The  nucleation  process  is  often  controlled 
either  by  the  stress  (tensile  pressure)  or  by  the  strain  state  experienced  by 
the  material.  For  instance,  in  the  biaxial  stretching  of  a  thin  sheet,  the 
process  is  controlled  by  plastic  deformation.  However,  in  a  plane  plate  im¬ 
pact  (uniaxial  strain)  condition,  the  void  nucleation  is  controlled  by  the 
high  tensile  stresses  generated  due  to  shock  wave  interactions.  The  nucleation 
threshold  strain  under  a  uniaxial  stress  state  is  always  orders  of  magnitude 
higher  than  under  uniaxial  strain.  In  general,  the  high  triaxial  stress  state 
enhances  the  ductile  void  growth  process.  The  modeling  of  the  ductile  failure 
process  requires  careful  consideration  of  the  influence  and  sensitivity  of 
stress/strain  states  as  well  as  strain  rates  on  the  nucleation  and  growth 
processes  of  microvoids. 

The  main  objective  of  this  report  is  to  describe  the  recently  developed 
dynamic  failure  model  [1]  and  its  applications  in  detail.  However,  prior  to 
describing  this  model,  the  authors  provided  a  comprehensive  background  on 


dynamic  failure  models  for  ductile  metals  in  Section  2.  The  essential  equa¬ 
tions  of  the  model  are  given  in  Section  3.  A  few  application  problems  are 
described  in  Section  4.  For  example,  spallation  in  a  target  due  to  a  long  rod 
penetration,  necking  in  a  dynamically  stretched  tensile  specimen,  a  spalla¬ 
tion  in  a  solid  cone  impacted  by  a  flyer  plate  on  its  base  are  simulated  using 
the  failure  model.  In  the  summary  section,  the  salient  features  of  the  new 
failure  model  are  discussed  and  recommendations  for  additional  experiments  and 
computer  simulations.  The  numerical  solution  scheme  of  the  governing  equa¬ 
tions  is  discussed  in  Appendix  A  and  B. 
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2.  BACKGROUND 


Numerous  journal  articles,  technical  reports,  and  conference  proceed¬ 
ings  on  dynamic  failure  models  exist.  Investigators  with  widely  varying 
backgrounds,  such  as  applied  mechanicists ,  metallurgists,  and  shock 
physicists,  have  contributed  to  the  understanding  of  dynamic  ductile  failure. 
Based  on  the  metallurgical  evidence,  Argon  et  al.  [2]  suggested  that  one  of 
che  most  often  observed  void  nucleation  processes  is  the  debonding  of  the 
matrix  material  from  the  second  phase  partictes.  Several  other  metallurgical 
investigations  support  this  evidence  [3-5],  McClintock  [6]  proposed  theories 
describing  growth  of  cylindrical  voids  and  their  subsequent  coalescence.  Rice 
and  Tracey  [7]  considered  the  growth  kinetics  of  a  single  spherical  void  in  an 
infinite  matrix  material  and  described  mathematically  the  dependence  of  void 
growth  on  triaxiality  of  the  stress  state.  Gurson  [8]  derived  a  pressure  de¬ 
pendent  yield  condition  for  porous  metals.  Using  this  condition,  he  proposed 
a  void  growth  law  similar  to  Rice  and  Tracey's.  Hancock  and  MacKenzie  ( 9  j 
provided  che  experimental  evidence  to  support  this  dependency  on  the  stress 
state  and  used  this  model  to  predict  failure  initiation  in  notched  tensile 
specimens  under  quasi -static  loading  conditions.  However,  the  mean  stress 
(pressure)  level  in  their  study  is  relatively  low  when  compared  with  the  le.-el 
in  a  plate  impact  test  configuration.  Also,  the  strain  rate  effects  become 
important  under  impact  loading  conditions. 

In  a  plate  impact  test,  tension  arising  from  the  interactions  of 
reflected  shock  waves  from  the  stress  free  planes  parallel  to  the  plane  of  im¬ 
pact  induces  fracture  in  the  target.  Void  nucleation  and  growth  occur  under 
very  high  triaxial  stress  leading  to  the  generation  of  a  stress  free  plane  in¬ 
side  the  target  through  coalescence  of  the  voids.  This  type  of  failure  is 
often  referred  to  as  spallation.  A  simple  criterion  in  which  stress  or  pres¬ 
sure  is  assumed  to  reach  a  critical  value  is  usually  employed  to  model 
spallation  in  computer  code  calculations.  This  type  of  time  independent 
failure  model  may  work  well  when  the  spall  is  well  abovj  the  threshold  condi¬ 
tions.  Tuler  and  Butcher  [10]  proposed  a  time  dependent  spall  criterion  in 
which  the  integrated  value  of  a  stress  dependent  damage  parameter  is  assumed 
to  reach  a  critical  value.  The  criterion  assumes  that  spall  occurs  if  the 
material  is  subjected  to  a  stress  greater  than  some  threshold  level  for  a  suf¬ 
ficient  characteristic  time.  In  this  model,  the  damage  parameter  does  not 
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degrade  either  the  stiffness  or  the  strength.  This  model  has  three  parameters 
which  have  to  be  determined  from  the  spall  experiments.  Rajendran  and  Bless 
[11]  determined  the  Tuler- Butcher  model  constants  for  several  metals. 

Cochran  and  Banner  [12]  proposed  a  relatively  advanced  model  to 
describe  the  spallation.  This  model  incorporates  the  effects  of  void  nuclea- 
tion  and  growth  on  the  strength  calculation  in  a  simplistic  manner.  The 
stress  is  allowed  to  relax  due  to  damage  evolution.  The  main  shortcoming  of 
this  model  is  that  it  is  limited  to  only  rectilinea,  motions  (one  dimensional 
strain).  Steinberg  et  al .  [13]  successfully  applied  this  model  in  conjunction 
wich  viscoplastic  constitutive  models  to  describe  the  spall  phenomenon  in 
several  metals.  The  dynamic  failure  process  in  metals  was  extensively  modeled 
using  a  microstatistical  approach  by  Curran  et  al.  [14-16].  This  relatively 
comple.;,  but  conceptually  advanced,  computational  model  pioneered  the  statis¬ 
tical  nature  of  the  nucleation  process  and  the  pressure  driven  void  growth 
process.  They  also  phenomenologically  modeled  the  coalescence  process.  This 
model  contained  a  number  of  parameters  determined  from  the  direct  measurements 
of  voids  in  the  target  of  a  carefully  conducted  recovery  experiment. 

Therefore,  the  determination  of  the  failure  model  parameters  is  unfortunately 
complex  and  time  consuming.  The  constitutive  equation  of  the  matrix  material 
is  relatively  simple  and  the  strength  degradation  is  handled  in  an  ad  hoc  man¬ 
ner,  rather  than  through  a  continuum  mechanics  based  three  dimensional 
approach.  However,  the  effects  of  damage  on  stiffness  and  strength  are  incor¬ 
porated  in  a  manner  which  leads  to  numerically  efficient  solutions. 

Recently  Curran,  Seaman,  and  Shockey  [16]  described  their  model  in 
detail  and  provided  an  in-depth  investigation  of  the  dynamic  failure  process 
in  both  ductile  and  brittle  metals.  Davison  et  al .  [17]  presented  a  theorv 
of  spall  damage  based  on  a  unified  and  thermodynamically  consistent  treatment 
of  elastic -viscoplastic  deformation.  A  continuum  mechanics  approach  was 
employed  with  a  three  dimensional  theory.  They  considered  a  detailed  disloc.; 
tion  based  thermoviscoplastic  constitutive  model  for  the  matrix  material 
behavior.  This  complete  treatment  jf  material  flow  and  failure  process  led  • 
numerous  model  constants  which  are  not  easily  obtainable. 

Johnson  and  Addessio  [18]  described  their  approach  to  model  tensile 
plasticity  and  ductile  fracture.  They  proposed  a  model  valid  for  general  ter 
sile  loading  conditions.  Therefore,  this  model  focused  on  describing  plastic 
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flow  and  fracture  in  notched  and  smooth  bars  as  well  as  spallation  under 
uniaxial  strain  conditions.  For  a  simplistic  approach,  the  nucleation  of 
voids  is  assumed  to  occur  instantaneously,  and  a  small  initial  void  content  is 
assigned  to  the  void  volume  fraction.  The  plastic  flow  was  described  by  the 
associative  flow  rule  derived  from  Gurson's  yield  surface  [8],  The  growth  law 
reduces  to  an  equation  derived  by  Carrol  and  Holt  [19]  for  tensile  threshold 
pressure.  The  matrix  material  behavior  is  described  through  a  rate  independ¬ 
ent  perfectly  plastic  model.  The  plastic  flow  in  the  porous  aggregate  is  also 
rate  independent.  Numerical  results  from  three  different  (one  dimensional 
strain,  two  dimensional  stress,  and  one  dimensional  stress)  calculations  were 
presented  to  demonstrate  the  model  capabilities. 

Rajendran  [20],  and  Rajendran,  Dietenberger ,  and  Grove  [1]  presented  a 
three  dimensional,  continuum  mechanics  based  dynamic  failure  model  (RDG  Model) 
to  describe  spallation  in  two  dimensional  stress  (long  rod  impact  on  a  thick 
plate)  and  one  dimensional  strain  (plane  plate  impact).  The  RDG  model  con¬ 
sidered  a  viscoplastic  constitutive  description  for  the  matrix  and  the  porous 
aggregate  materials.  The  stress  and  strain  based  void  nucleation  process  was 
modeled  through  a  Gaussian  distribution,  as  initially  proposed  by  Chu  and 
Needleman  [21]  to  describe  strain  rate  independent  nucleation  process  in  thin 
sheet  stretching.  The  RDG  model  proposed  a  new  pressure  dependent  yield  func¬ 
tion  for  describing  plastic  flow  in  the  porous  aggregate.  There  are  four 
phases  in  the  model.  In  the  first  phase,  the  intact  material  is  described  by 
the  Bodner- Partom  viscoplastic  model  [22].  The  void  nucleation  is  introduced 
in  the  second  phase.  The  void-contained  aggregate  is  described  in  the  third 
phase  using  an  associated  plastic  flow  rule  derived  from  a  pressure  dependent 
yield  function.  The  last  phase  of  modeling  is  the  coalescence  of  voids  lead¬ 
ing  to  complete  failure.  In  the  RDG  model,  separate  modeling  of  the 
coalescence  process  is  not  needed.  The  void  growth  law  is  such  that  the 
growth  rate  is  rapidly  increased  as  the  damage  approaches  its  critical  value. 
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3.  FAILURE  MODEL  FORMULATION’ 


Three  of  che  main  features  in  a  dynamic  ductile  failure  model  are  the 
following:  one,  the  initial  intact  (void-free)  material  requires  a  constitu¬ 

tive  description  which  will  include  strain  rate  and  temperature  effects;  two, 
a  mathematical  description  of  the  void  nucleation  and  growth  process;  and 
three,  pressure  dependent  plastic  flow  equations  for  the  porous  (void- 
containing)  aggregate.  There  could  be  several  other  features  in  the  failure 
modeling.  For  instance,  depending  on  the  level  of  voids,  the  stiffness  of  the 
aggregate  may  significantly  be  reduced;  therefore,  models  to  degrade  the 
stiffness  will  be  required  for  the  realistic  description  of  the  material  be¬ 
havior.  The  final  process  of  voids  coalescence  may  also  be  important. 

Another  aspect  of  the  failure  model  is  the  equation  of  state  for  the  porous 
aggregate. 

A  brief  description  of  the  Bodner-Partom  (BP)  model  for  the  intact 
(matrix)  material  is  given  in  this  section.  As  an  alternative  to  BP  model, 
the  model  proposed  by  Johnson  and  Cook  (JC)  [23]  is  also  described.  The  model 
constants  for  BP  and  JC  are  provided  for  several  metals.  The  model  constant 
evaluation  scheme  for  these  two  models  have  been  described  elsewhere, 
reference  [24]  for  the  BP  model  and  reference  [25]  for  JC  model. 

3.1  Bodner-Partom  Model 

The  fu' ly  densed,  void- free  matrix  material  can  be  modeled  through  the 
sLdle  variable  based  viscoplastic  constitutive  equations  of  Bodner  and  Fartom 
(B-P  modal)  [22].  The  B-P  model  in  terms  of  equivalent  plastic  strain  rate, 

•  p 

D  ,  and  effective  stress,  i  ,  is  given  by, 
ra  m  ° 


_2. 

73 


m 


2n 


(1) 


where  Z  is  a  state  variable.  D  is  the  limiting  value  of  the  plastic  strain 

°  8 

rate.  The  value  of  Dq  is  usually  set  to  10  /sec  for  metals.  n  is  a  parameter 
that  is  mainly  related  to  strain  rate  sensitivity.  The  state  variable  Z 
describes  the  overall  resistance  of  the  material  to  plastic  flow,  and  it 
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depends  on  che  loading  history.  The  evolution  equation  for  the  loading  his¬ 
tory  parameter  Z  is 

Z  -  m(Z,  -  Z)  V  (2) 

1  P 

• 

where  W  is  the  plastic  work  rate.  Z.  is  the  maximum  value  that  Z  can  attain, 
P  1 

and  m  is  a  parameter  that  embodies  the  strain  hardening  behavior  of  the 
material.  For  highly  strain  hardening  materials  like  copper,  m  is  described 

by, 


-aW 

m  -  m  +  m.,  e  ^  (3) 

o  1 

where  m^,  m^  and  a  are  additional  model  parameters.  For  less  strain  hardening 
materials,  m^  and  a  are  assumed  to  be  zero.  For  example,  the  modeling  of  OFHC 
copper  required  the  additional  two  constants. 

The  high  temperature  effect  on  the  strength  (thermal  softening)  is  in¬ 
corporated  simplistically  through  the  strain  rate  sensitivity  constant,  n. 
Experimental  evidence  suggests  that  at  elevated  temperatures  metals  exhibit 
increased  rate  dependency.  The  temperature  dependent  n  is  described  by, 

n  -  A  +  B  /  T°  (4) 

where  the  temperature  is  expressed  in  degrees  Kelvin.  A  and  B  are  additional 

BP  model  constants.  Previously,  Rajendran  et  al .  [24]  presented  a  systematic 

scheme  for  evaluating  the  model  constants.  They  determined  the  B-P  model 

parameters  from  split  Hopkinson  bar  experimental  data  and  the  steady  state 

value  of  cr  .  The  BP  model  constants  for  several  metals  have  been  determined 
H&L 

and  tabulated  in  Table  1.  Z„  is  the  initial  value  of  Z.  The  D  is  assumed  to 

0  8  ° 

be  the  same  for  all  metals,  and  a  value  of  10  /sec  has  been  arbitrarily  as¬ 
signed.  In  his  review  report,  Bodner  [26]  suggested  a  value  of  10^/sec  based 
on  dislocation  motion  measurements  for  high  velocity  impact  loading  condi¬ 
tions.  A  lower  value  for  the  strain  rate  sensitivity  parameter  n  indicates 
that  the  material  is  extremely  strain  rate  sensitive.  The  C1008  steel  is  an 
example.  The  value  of  4  for  6061-T6  aluminum  indicates  that  the  material  is 
the  least  strain  rate  dependent.  When  the  difference  between  Z^  and  Z^  is 
small,  the  material  is  relatively  less  strain  hardening.  For  OFHC  copper  the 
difference  is  significantly  high  which  means  high  strain  hardening  behavior. 
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TABLE  1 

BODNER-PARTOM  MODEL  CONSTANTS 


Material 

Zo 

(GPa) 

z, 

(GPa) 

n 

m0 

GPa1 

m, 

GPa1 

a 

GPa1 

A 

B 

Cl  008 

Steel 

5.5 

7.0 

0.4 

15 

0 

0 

0.245 

46 

HY100 

Steel 

2.4 

3.6 

H 

10 

0 

0 

NA 

NA 

1020 

Steel 

0.64 

0.93 

4.0 

30 

0 

0 

NA 

NA 

MAR-200 

Steel 

2.2 

2.4 

4.0 

5 

0 

0 

NA 

NA 

Armco 

Iron 

2.65 

4.2 

0.58 

56 

0 

0 

NA 

NA 

OFHC 

Copper 

0.8 

6.6 

0.4 

11 

150 

1500 

NA 

NA 

6061-T6 

Aluminum 

0.45 

0.55 

4.0 

120 

0 

0 

-2.86 

2343 

7039-T64 

Aluminum 

0.56 

0.76 

4.0 

28 

0 

0 

NA 

NA 

Pure 

Tantalum 

— 

3.1 

0.74 

20 

0 

0 

NA 

NA 

W-2 

Tungsten 

8.75 

10.0 

0.58 

150 

0 

0 

0.166 

134 

Nickel 

200 

0.32 

0.82 

4.0 

40 

0 

0 

NA 

NA 

MAR-250 

Steel 

2.5 

2.7 

5.0 

20 

0 

0 

NA 

NA 

AF1410 

Steel 

■ 

2.75 

5.0 

15 

0 

0 

NA 

NA 

NA  --  The  high  temperature  constants  are  "Not  Available" 
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3.2  Johnson-Cook  Model 


The  Johnson-Cook  constitutive  model  was  used  as  a  basis  of  comparison 
for  this  study.  This  model  is  basically  an  empirical  model  which  incorporates 
strain  rate  and  temperature  dependency  on  strength.  The  plastic  strain  rates 
are  not  bounded  and  therefore  are  not  restricted  below  any  limiting  value. 
Also,  no  loading  history  effects  are  built  into  the  model.  However,  this 
five -parameter  empirical  model  has  been  demonstrated  to  provide  realistic 
solutions  to  hydrocode  predictions  of  a  very  broad  class  of  applications  to 
extremely  dynamic  events  such  as  impact,  penetration,  and  explosive  accelera¬ 
tion  of  metals.  The  success  of  the  model  is  due  to  its  ability  to  incorporate 
the  overall  effects  of  temperature  and  strain  rate  as  observed  from  the 
dynamic  and  quasi-static  experimental  data.  The  Johnson-Cook  model  has  the 
following  form: 


Y  -  [A  +  B  £n]  [l  +  C  in  7*]  [l  -  T*M]  (5) 


* 


T 


T  -  T 
_ room 


melt 


T 

room 


(6) 


where  Y  is  the  flow  strength,  «  is  effective  plastic  strain,  e*  is  non- 
dimensional  strain  rate  (normalized  by  1/sec) .  The  five  material  constants 
are  the  following:  A  is  the  yield  strength,  B  is  the  work  hardening  coeffi¬ 
cient,  n  is  the  work  hardening  exponent,  C  is  strain  rate  coefficient,  and  M 
is  the  thermal  softening  exponent.  In  Table  2,  the  values  of  these  constants 
are  tabulated  for  various  metals. 


The  model  constants  A,  B,  and  n  are  determined  from  quasi-static 
stress  -  strain  data  either  from  tensile  or  compressive  tests.  The  strain  rate 
dependent  constant  C  can  be  determined  from  the  slope  of  the  stress  versus 
strain  rate  plot.  The  temperature  constant  M  is  estimated  from  the  stress 
versus  temperature  plot.  In  hydrocode  analysis,  strain  rate,  and  the  von- 
Mises  yield  radius  can  be  calculated  from  the  JC  equation  for  a  given  strain 
and  temperature. 

3.3  Nucleation  and  Growth  of  Voids 


Until  voids  nucleate,  the  aggregate  behavior  can  be  described  by  the  F 
P  model.  The  plastic  flow  in  the  void- free  aggregate  is  incompressible,  i.e 
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STRENGTH  MODEL  CONSTANTS  FOR  JOHNSON-COOK  MODEL 
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not  provided 


the  sum  of  the  principal  or  orthogonal  plastic  strains  is  equal  to  zero. 
However,  the  nucleation  of  voids  will  introduce  dilatation  and  the  plastic 
yield  behavior  will  depend  not  only  on  the  second  invariant,  J2,  but  also  on 
the  mean  stress  or  pressure.  The  constitutive  model  for  the  aggregate  must 
include  these  effects.  For  this  purpose,  a  yield-criterion-based  plastic  flow 
rule  in  which  the  pressure  dependence  enters  explicitly  into  the  calculations 
was  considered. 

The  most  widely  used  void  nucleation  model  was  the  one  that  was 
initially  used  by  Chu  and  Needleman  [21]  in  their  analysis  of  localized  neck¬ 
ing  in  biaxinlly  stretched  sheets.  The  model  was  based  on  a  mechanism  in 
which  voids  are  nucleated  due  to  debonding  of  inclusions  from  the  matrix.  The 
debonding  can  occur  because  of  a  stress  or  a  strain  criterion  or  both.  The 
corresponding  model  for  void  nucleation  rate  is  given  by, 


where 


and 


f  -  F  (Y  +  P)  +  F  DP 
n  a  m  e  m 


F 

a 


h  In 

2[ 


P 

S1 


2 

) 


F 

€ 


(7) 


(8) 


(9) 


If  the  nucleation  is  due  to  cnly  the  matrix  debonding  from  inclusions, 
then  the  total  void  volume  nucleated  must  be  consistent  with  the  volume  frac- 
:ion  of  second  phase  particles.  Therefore,  the  values  determined  for  the 
parameters  f^  and  must  meet  this  requirement.  and  e^  are  the  mean 

equivalent  stress  and  strain,  respectively,  around  which  the  nucleation  stress 
and  strain  are  distributed  in  a  Gaussian  manner.  s^  and  S2  are  the  standard 
deviations  of  the  distributions.  These  two  parameters  will  control  the  ranges 
of  stress  or  strain  over  which  most  of  the  voids  can  be  nucleated. 


11 


The  growth  law  can  be  directly  related  to  the  dilatation  due  to  growth 
of  voids  in  tne  aggregate.  By  definition,  the  void  growth  rate  is  given  by, 


f 

g 


P 


* 


(10) 


•  p 

where  repeated  index  means  summation,  and  are  plastic  strain  rates  in  the 
three  principal  directions  ana  f  -  1-/3,  where  p  is  the  relative  density,  p  is 
the  ratio  of  current  density  to  the  reference  density. 


3.4  Pressure  Dependent  Plastic  Flow  of  the  Aggregate 


A  pressure  dependent  yield-criterion-based  approach  has  been  considered 
in  the  constitutive  model  formulation.  For  randomly  distributed  voids  or 
microcracks  contained  in  the  aggregate,  the  yield  behavior  will  be  influenced 
by  not  only  the  second  invariant  of  the  deviatoric  stress  (J2)  t>uC  a^so  by  che 
pressure  or  mean  stress  (1^).  The  following  form  of  the  yield  function  has 
been  considered: 


A (p>  J2  +  B(p)  I*  -  S(p)  Y*  (11) 

where  A,  B,  and  S  are  functions  of  relative  density,  p.  Y^  is  the  effective 
stress  in  the  matrix  material.  (Note:  the  subscript  'm'  means  matrix  material 
and  not  a  Censorial  index.)  Based  on  a  critical  total  deformation  energy, 
Doraivelu  et  al.  [27]  derived  the  following  expressions  for  A  and  B: 


A (p)  >  2  +  p2 

(12) 

2 

B(P)  -  L  3  * 

(13) 

In  general,  the  5  function  is  material  dependent  while  the  functions  A  and  B 
are  independent  of  the  matrix  material  behavior.  The  conditions  on  the  coef¬ 
ficient  of  the  matrix  effective  stress  in  the  yield  function  are  5(p)  -  1  at 
-  1  and  S(p)  -  0  at  p  -  Pcx-  Rajendran  et  al.  [1]  proceeded  initially  with 
the  form 


g (P)  -  g (Prr) 

‘  «(i>  •  *<,„> 


(14) 


12 


which  satisfies  the  conditions  on  the  6  function.  Rajendran  et  al.  also 
proposed  a  new  function  for  g (p)  with  the  form 


S (P)  "  [<«  *  ]f[  (i  -  P)  f  (15) 

where  k  and  N  are  model  constants.  A  negative  value  of  N  makes  Equation  (14) 

a  hyperbolic  power  function.  Numerical  simulations  of  a  plate  impact  test 

configuration  were  performed  to  evaluate  the  effects  of  the  5  function  on  the 

space  signal.  In  these  simulations,  the  initial  slope  of  S(p)  at  p-1  great  iy 

influenced  the  slope  of  the  spall  signal.  By  taking  advantage  of  this  aspect, 

a  parameter,  /?  (-S'(l))  was  introduced.  With  /9,  N,  and  p  as  model  con- 

cr 

stants,  the  corresponding  value  of  k  can  be  solved  by  a  simple  iterative 

scheme.  An  idealistic  value  of  one  can  be  assumed  for  o  Therefore.  N  and 

cr 

/3  become  the  only  model  parameters  for  void  growth  description. 

As  a  first  step  in  deriving  the  plastic  strain  rates,  the  yield  condi¬ 
tion  for  the  aggregate  can  be  rewritten  as: 

2 

*  -  (2+p2)  J2  +  l\  ■  6(p)  Y2  -  0  (16) 

The  viscoplastic  strain  rates  in  the  aggregate  can  be  calculated  using 
the  flow  rule  derived  as: 


•  3$ 

do.  . 
1J 


(17) 


The  proportionality  factor,  A,  can  be  obtained  using  the  flow  rule  in  conjunc¬ 
tion  with  the  following  relationship: 


(1-f)  Y  DP 
m  m 


(18) 


where  f  is  the  void  volume  fraction  and  related  to  relative  density  p  through 
f  -  1 -p.  Note  that  'm'  is  a  subscript  and  not  a  tensorial  index.  The  above 
expression  was  derived  from  the  definition  that  the  total  plastic  work  in  the 
aggregate  is  entirely  due  to  the  plastic  work  done  by  the  matrix.  By  combin- 
ing  Equations  (17)  and  (18) ,  A  can  be  expressed  by, 
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(1-f) 


X  - 


Y 

m  m 


a$ 

da.  . 

1J 


ij 


(19) 


The  plastic  strain  rates  in  the  aggregate  can  be  rewritten  as: 


(1-f)  Y  DP 
_ m  m 

da  arl 
rl 


(20) 


In  the  constitutive  model  formulation,  the  total  strain  rate  is  decomposed 

into  elastic  and  plastic  strain  rates.  The  elastic  strain  rates,  «...  are  re- 

ij 

lated  to  the  stress  through  Hooke's  law  as: 


Dik 


(21) 


where  D.,  is  the  inverse  of  elastic  modulus  matrix,  E.,  . 

lk  lk 

Using  the  consistency  condition  which  holds  during  the  plastic  flow,  an 
expression  can  be  obtained  for  Y^  as: 


Y 

m 


r  •  aj> 
^  da . .  aij  +  df 
— U - - 

dY 


n 


(22) 


An  expression  for  a .  .  can  be  obtained  from  Equation  (21)  by  replacing  the 
elastic  strain  rate  as  the  difference  between  total  and  plastic  strain  rates. 


a .  . 


(23) 


3.5  Degradation  of  Shear  and  Bulk  Moduli 

The  shear  and  bulk  moduli  are  degraded  using  Mackenzie's  [28]  formula 
Mackenzie  derived  an  approximate  analytical  expression  for  the  reduction  of 
elastic  stiffness  based  on  the  elastic-plastic  flow  around  a  spherical  void 
an  infinite,  incompressible  matrix  material.  The  RDG  model  [1]  employed  the 

A  A 

corresponding  expression  for  the  degraded  moduli  K  and  G  as, 
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a  -  n 


(i  + 


3Kf 

4G 


(24) 


G_  (1  -  f)  (1  -  £6K  ±  12G)  f)  ,  (25) 

G  (9K  +  8G) 


where  K  and  G  are  the  bulk  and  shear  moduli  of  the  intact  material.  Since 
the  void  volume  fraction  is  also  changing  with  time,  the  stress  rate  equation 
(23)  needs  to  be  modified  to  add  a  term  proportional  to  the  void  growth  rate. 
Appendices  A  and  B  outline  the  solution  scheme  to  solve  the  governing  equa¬ 
tions  . 


3.6  Equation  of  State  for  the  Aggregate 

The  equation  of  state  for  the  solid  (intact)  material  has  been 
routinely  determined  from  the  experimentally  obtained  Hugoniot  data  for 
several  metals,  such  as  steel,  aluminum,  copper,  etc.  However,  once  the 
initially  intact  solid  develops  voids  upon  tensile  loading,  the  equation  of 
state  (EOS)  of  the  original  solid  material  is  no  longer  valid  during  the  sub¬ 
sequent  compressive  loading.  For  this  purpose,  the  EOS  parameters  have  to  be 
modified. 

The  simple  approach  used  in  the  RDG  model  is  an  extension  of 
Mack  nzie's  procedure.  The  EOS  for  the  intact  solid  is  described  by  the  Mie- 
Gruneisen  equation  as 

Pc  -  (p  n  +  09m2  +  £,m3) d-rM/2)  +  r( i - i  )  ,  (26) 

S  1  2  i  O 


where  0^,  0^,  0^,  and  T  are  EOS  constants.  ^  is  the  volumetric  compressible 
strain  of  the  intact  solid,  (V  /V-  1)  ,  and  I  is  the  specific  internal  energy 
of  the  solid  material.  To  account  for  the  voids,  the  Mackenzie's  formulation 
is  implemented  as  follows.  A  volumetric  elastic  strain  of  the  aggregate  is 
first  defined  as, 


^ag-  [  (l+^/d-f)  -  1  ]  ,  (27) 

to  replace  n  in  the  right  side  of  Equation  (26).  Then  the  right  side  expres- 

A 

sion  of  Equation  (26)  is  multiplied  by  the  Mackenzie  correction  term,  K/K  ,  to 
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obtain  the  aggregate  pressure,  P.  Thus  a  voided  material  undergoing  compres¬ 
sion  will  load  along  a  degraded  bulk  modulus  at  the  aggregate  compressible 
strain.  At  high  enough  stress  levels,  the  void  will  collapse  (f-0),  and  thus 
the  modified  EOS  will  reduce  to  Equation  (26). 

The  various  governing  equations  can  be  rearranged  and  numerically  in¬ 
tegrated.  The  RDG  model  along  with  BP  or  JC  model  was  implemented  into  the 
EPIC2  finite  element  code  [29],  Special  purpose  subroutines  have  been 
developed  successfully .  The  corresponding  numerical  scheme  is  given  in 
Appendix  A. 
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4.  FAILURE  ANALYSES  USING  THE  RDG  MODEL 


The  strength  model  (BP  or  JC)  constants  are  determined  using  the  stress 
-  strain  data  from  split  Hopkinson  bar  (SHB)  tests.  These  data  may  comprise 
data  from  quasi-static  tensile,  SHB  tensile,  compressive,  and  torsional  tests. 
Rajendran  et  al.  [24]  described  a  combined  experimental  and  numerical  scheme 
for  evaluating  BP  model  constants.  They  also  included  the  HEL  (Hugoniot 
Elastic  Limit)  data  in  the  scheme.  Johnson  and  Holmquist  [30]  provided  a 
methodology  to  obtain  JC  model  constants.  The  BP  and  JC  model  constants  for 
various  metals  are  given  in  Tables  1  and  2. 

The  determination  of  the  failure  model  constants  requires  plate  impact 
test  data  in  terms  of  either  velocity  history  or  stress  history  at  the  back  of 
the  target.  Rajendran  and  Bless  [11,31]  reported  the  test  and  diagnostic 
techniques  as  well  as  data  for  several  metals.  Rajendran  et  al .  [32]  deter¬ 

mined  the  failure  model  constants  using  the  velocity  history  obtained  from  a 
plate  impact  test  on  OFHC  copper. 

Numerical  exercises  were  conducted  using  the  new  subroutines.  The  ex¬ 
ercises  were  based  on  a  plate  impact  test  simulation. 

4.1  Plate  Impact  Simulations 

Detailed  discussions  on  planar  plate  impact  test  technique  can  be  found 
in  References  33  and  34.  Determination  of  the  failure  model  parameters  is 
aided  through  the  simulation  of  plate  impact  tests.  Simulations  are  carried 
out  using  the  STEALTH  one-dimensional  finite  difference  code  [35],  The  effect 
of  each  model  parameter  on  the  spall  behavior  is  evaluated  from  the  simulation 
results.  A  methodology  to  determine  the  parameters  is  outlined.  The  model 
also  has  been  successfully  used  to  describe  spallation  in  OFHC  copper. 

4.1.1  Physical  Features 

Impact  of  a  2  mm  thick  copper  flyer  against  a  9  mm  thick  OFHC  copper 
target  is  modeled  using  the  STEALTH  code.  The  constitutive  and  failure  model.- 
described  in  Section  3  were  used  to  characterize  the  high  strain  rate  behavi  n 
of  copper.  Our  first  objective  was  to  demonstrate  the  experimentally  observ.  : 
important  physical  features  of  the  free  surface  velocity  profile.  In  Figure 
1,  results  from  two  different  simulations  are  shown.  In  the  first  simulation 
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important  physical  features  of  the  free  surface  velocity  profile.  In  Figure 
1,  results  from  two  different  simulations  are  shown.  In  the  first  simulation 
the  model  parameters  were  chosen  so  that  the  cargat  spalls.  A  typical  spall 
signal,  as  is  usually  observed  in  an  experiment,  can  be  seen  in  Figure  1.  In 
the  second  simulation  the  spall  is  suppressed  by  choosing  a  zero  value  for  the 
nucleation  parameter,  f^.  The  complete  unloading  of  the  velocity  history,  as 
indicated  by  the  dotted  line,  clearly  demonstrates  the  absence  of  spall  in  the 
target . 

In  Figure  2,  the  velocity  histories  obtained  from  simulations  at  dif¬ 
ferent  impact  velocities  are  shown.  The  most  important  physical  features  in 
th  •.  velocity  profiles  are  the  velocity  level,  AV^ ,  and  the  time  duration  or 
period,  r  ,  of  the  wave  transit  between  the  spai.  plane  and  the  free  surface. 
The  AV^  corresponds  to  a  stress  level  around  which  rapid  microvoid  or 

microcrack  nucleation  occurs  (note  that  a.,  -  1/2  p  C  AV  ) .  If  the  impact 

Ns 

velocity  is  greater  than  AV^ ,  as  in  the  cases  of  plots  A,  B,  C,  and  D  in 

Figure  2,  then  spallation  will  occur,  as  indicated  by  the  spall  signal. 

uowever,  at  an  impact  velocity  oi  50  m/s,  spall  nucleation  has  not  occurred. 

The  tensile  stresses  in  the  target  at  this  impact  velocity  are  lower  than  the 

mean  nucleation  threshold  stress  cr,  (=16  kbars) .  In  the  nucleation  model, 

N 

nucleation  is  assumed  to  occur  at  <7„  ±  3s,  where  s  is  the  standard  deviation. 

N 

Figure  3  shows  that  at  velocity  V  -  100  m/s,  the  tensile  peak  stress  was 
around  13  kbars,  compared  to  a  peak  stress  of  18  kbars  in  compression.  This 
is  due  to  void  softening  of  material  which  reduced  the  stress  levels. 

The  void  volume  fraction  distribution  levels  in  the  target  at  three 
different  times  are  shown  in  Figure  A.  The  distribution  clearly  shows  that 
the  maximum  void  volume  fraction  is  at  the  spall  plane  whi ~h  is  around  2  mm 
(flyer  plate  thickness)  from  the  free  surface  (at  x  -  7  mm).  Presence  of 
voids  at  and  around  the  spall  plane  has  been  supported  by  metallographical 
studies  conducted  on  different  materials  [14-16]. 

In  Figure  5,  the  loading  path  at  the  spall  plane  of  the  target  is 
shown.  The  void  volume  levels  are  shown  by  dotted  lines.  Initially,  the 
strength  is  independent  of  pressure  as  can  be  seen  between  points  A-B.  Damage 
nucleation  has  not  yet  initiated,  and  therefore  f  remairs  zero.  At  B,  the 
nucleation  occurs.  As  the  tensile  pressure  increases,  the  void  volume  also 
increases  between  the  loading  points  B-C.  Strength  rapidly  decreases  between 
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Figure  1.  Simulated  Velocity  History  for  a  Target  without  and  with  Spall  Failure 
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Figure  3.  Simulated  Stress  History  at  the  Spall  Plane.  Stress  Relaxation  under 
Tension  is  Shown  for  Three  Impact  Velocities  -  100,  185,  and  500  m/s. 
(Note:  stress  is  -Ve  in  tension  and  +Ve  in  compression.) 


:ion  in  the  Impacted 


points  C-D.  At  one  percent  void  (around  D) ,  the  void  contained  aggregate  can 
no  longer  sustain  tensile  pressure,  so  the  pressure  rapidly  decreases  as  the 
void  volume  reaches  10  percent.  Failure  (coalescence  of  voids)  occurs  between 
points  E-F.  Both  pressure  and  strength  approach  zero  as  the  material  com¬ 
pletely  fails. 

4.1.2  Stability  of  the  Solution 

The  integration  of  the  ordinary  differential  equations  for  the  con¬ 
stitutive  model  by  a  Diagonally  Implicit  Runge-Kutta  (DIRK)  scheme  is  known  to 
be  stable  and  at  least  second  order  accurate.  But  if  the  time  intervals  of 
the  DIRK  scheme  are  much  smaller  than  the  time  step  determined  by  the  Courant 
criterion  in  the  STEALTH  wave  propagation  code  [35],  then  the  stability  and 
accuracy  of  the  solution  becomes  uncertain.  Depending  on  the  DIRK  time  step 
size,  the  STEALTH  stable  time  step  for  a  particular  zone  (or  element),  as 
determined  from  the  Courant  condition,  may  be  reduced  by  an  integer  factor 
ranging  from  2  to  10.  This  factor  is  initialized  to  one,  indicating  no  reduc¬ 
tion  in  the  stable  time  step.  Then  the  factor  is  either  increased  by  one,  if 
the  last  DIRK  time  step  is  less  than  one  tenth  of  the  STEALTH  time  step,  or 
decreased  by  one,  if  the  last  DIRK  time  step  is  greater  than  two  tenths  of  the 
STEALTH  time  step.  This  method  allows  a  gradual,  but  not  necessarily  per¬ 
manent,  reduction  in  the  STEALTH  stable  time  step  down  to  one  tenth  of  the 
stable  time  step  computed  from  the  Courant  condition.  Since  there  is  no  math¬ 
ematical  criterion  of  stability  for  the  entire  solution,  a  few  numerical 
exercises  were  devised  to  check  for  solution  stability  and  accuracy.  These 
exercises  involved  several  plate  impact  simulations.  One  effective  numerical 
test  was  to  determine  if  the  total  momentum  and  the  total  energy  remained  con¬ 
stant  during  a  plate  impact  simulation.  The  solution  was  accepted  only  if 
this  condition  was  met.  The  second  exercise  was  to  determine  if  the  entire 
solution  remained  essentially  the  same  if  the  grid  size  was  varied.  The 
results  are  shown  in  Figure  6  for  two  impact  velocities  (200  m/s  and  500  m/s) 
with  three  different  grid  sizes.  It  is  possible  to  ensure  stability  and  ac¬ 
curacy  by  forcing  the  STEALTH  time  interval  to  be  the  same  as  the  DIRK  scheme 
time  interval,  but  this  is  costly  in  terms  of  programming  effort  and  excessive 
computer  time.  In  any  case,  numerical  exercises  to  vary  the  numerical  in¬ 
tegration  and  stability  parameters  should  be  conducted  to  optimize  the 
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Fiqurt:  5.  The  Loading  Path  (/3J2/Ym  vs.  3P/Ym)  at  the  Spall  Plane  in 


stability,  accuracy,  and  computer  time  of  the  solution.  The  results  so  far 
indicate  that  the  solution  is  stable  and  reasonably  accurate. 

4.1.3  Sensitivity  Study 

It  is  important  to  evaluate  the  effects  of  various  model  parameters 
(a^,  f^,  s,  p ,  N,  f£r)  on  the  numerically  simulated  failure  process.  For  this 
purpose,  a  sensitivity  study  was  considered.  The  sensitivity  of  the  dynamic 
failure  model  parameters  on  the  solution  (in  the  free  surface  velocity  versus 
time  plot)  can  be  checked  by  varying  the  values  systematically.  The  variation 
in  the  values  of  the  three  void  nucleation  parameters,  ( a ,  f^  and  s)  were  ex¬ 
amined  first,  and  then  the  three  void  growth  parameters  (p,  N,  and  f  )  were 
examined. 

An  increase  in  the  value  of  caused  the  spall  signal  to  occur  later 

while  the  solutions  during  the  rebound  tend  to  merge  together.  Likewise,  a 

decrease  in  the  value  of  f^  caused  the  spall  signal  to  occur  later  while  the 

solutions  during  the  rebound  tend  to  merge  together.  It  appears  the  void 

nucleation  parameters  and  f^  have  a  strong  influence  on  the  initial  spall 

signal,  but  become  a  negligible  influence  on  the  void  growth  which  affects  the 

rebound.  The  standard  deviation  s  can  be  defined  as  a  fraction  of  c\T.  ”*he 

N 

sensitivity  study  used  three  different  fractions,  0.125,  0.25,  and  0.5. 

Recall  from  Equations  (7)  and  (8)  that  the  model  for  void  nucleation  due  to 

stress  follows  a  Gaussian  distribution  with  mean  a.,,  which  allows  nucleation 

N 

to  occur  for  a  stress  range  of  cr^  ±  3s.  Since  void  nucleation  in  metals  can 
only  occur  during  tension,  a  practical  upper  limit  for  the  fraction  is  one 
third.  This  prevents  a  negative  lower  bound  for  the  nucleation  stress  range. 
Varying  the  stress  standard  deviation,  s,  had  only  a  minor  effect  on  the  spall 
signal  or  the  rebound. 

The  rebound  peak  of  the  spall  signal  increased  for  increasing  values  o: 
the  void  growth  parameters  p  and  |N|.  The  initial  slope  of  the  spall  signal 
was  less  influenced  by  N  than  p.  Further  sensitivity  studies  showed  that  the 
void  growth  parameter  P  seems*  to  affect  the  slope  significantly.  Taking  ad¬ 
vantage  of  the  negligible  influence  of  N  on  the  slope,  the  parameter  P  can  be 
estimated  by  matching  the  spall  signal  slopes  between  the  simulation  and  ex¬ 
periment.  Varying  f  from  0.6  to  1.0  had  a  very  minimal  impact  on  the  spall 
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Figure  6.  Demonstration  of  Grid  (or  mesh)  Independency  of  the  Solution  for 
Three  Different  Grids  at  Two  Velocities. 


In  summary,  the  sensitivity  study  indicates  that  the  failure  model 
parameters  can  be  systematically  determined  by  matching  the  simulation  results 
with  the  spall  signal.  Preliminary  estimate  for  can  be  obtained  from  the 
relationship,  a -  1/2  p  C  AV^ ,  where  p  is  the  mass  density  and  C  is  the  wave 
speed.  The  value  for  AV^  is  available  from  the  velocity  history.  The  value 
for  s  can  be  arbitrarily  chosen  to  be  one-fourth  of  a^.  A  theoretical  value 
of  1  can  be  assigned  for  f  .  The  only  remaining  parameters  that  need  to  be 
determined  from  the  spall  signal  are  f  ^ ,  £,  and  N.  f^  can  be  determined  by 
matching  the  spall  arrival  time;  0  is  chosen  based  on  matching  the  slope  of 
the  spall  signal;  finally,  N  is  selected  based  on  matching  the  rebound  peak  of 
the  spall  signal.  Rajendran  et  al.  [32]  employed  this  scheme  and  determined 
the  model  parameters  for  OFHC  copper.  The  following  section  describes  this 
effort . 

4.2  Model  Parameters 

For  establishing  a  standard  procedure  to  calibrate  the  failure  model 
constants,  two  plate  impact  experiments  (#538  and  #1299)  were  considered.  The 
first  experiment  was  on  an  OFHC  copper  target  plate  and  the  velocity  history 
was  obtained  using  VISAR  (velocity  interferometry) .  The  second  experiment 
considered  an  HY100  steel  target,  and  in  this  experiment  the  diagnostic  was 
the  manganin  gage  measured  stress  history.  The  test  configuration  details  are 
given  in  Table  3.  These  two  experiments  were  simulated  using  the  EPIC  code. 
The  special  purpose  subroutines  describing  the  RDG  model  were  used  to  model 
the  spall  process  in  both  OFHC  copper  and  HY100  steel.  Using  the  suggested 
procedure,  the  model  constants  for  these  two  as  well  as  several  other  metals 
were  determined.  The  corresponding  constants  are  given  in  Table  4. 

The  BP  model  was  employed  for  the  intact  material  description.  The  cor¬ 
responding  parameters  were  given  earlier  in  Table  1.  The  plate  impact 
experiment  on  annealed  OFHC  copper  reported  by  Rajendran  et  al.  [32]  was 
employed  in  the  failure  model  parameter  evaluation.  A  2  mm  copper  flyer  was 
impacted  against  a  9  mm  copper  target  at  an  impact  velocity  of  185  m/s.  The 
free  surface  velocity  history  of  the  target  was  determined  by  VISAR  measure¬ 
ments.  This  velocity  history  was  used  in  the  model  parameter  calibrations. 

Based  on  the  sensitivity  study  discussed  in  the  preceding  section,  the 
failure  model  parameters  are  estimated.  Out  of  the  six  parameters  (aN,  s, 


TABLE  3 

PLATE  IMPACT  EXPERIMENTS 


Shot 

Number 

Flyer  Plate 

Target  Plate 

Impact 

Velocity 

(m/sec) 

Material 

Thickness 

(mm) 

Material 

Thickness 

(mm) 

7-1299 

Copper 

3 

HY100  Steel 

6 

489 

7-1300 

Copper 

3 

C1008  Steel 

6 

462 

7-1298 

Copper 

3 

Armco  Iron 

6 

470 

7-1268 

1020  Steel 

3.9 

1020  Steel 

7.8 

572 

7-1288 

Copper 

2 

Copper 

4.1 

512 

7-1302 

Copper 

2 

Tantalum 

6 

735 

7-1455 

MAR-200  Steel 

3 

MAR-200  Steel 

6 

508 

7-1454 

MAR-250  Steel 

3 

MAR-250  Steel 

6 

575 

7-1523 

AF1410  Steel 

4.5 

AF1410  Steel 

9.0 

678 

TABLE  4 

RDG  MODEL  CONSTANTS 


Material 

Nucleation 

Yield  Function 

ft 

<*N 

s=(oN/4) 

P 

N 

OFHC  Copper 

0.01 

1.6 

0.4 

65 

-2.4 

HY100  Steel 

0.035 

3.6 

0.9 

30 

-1.0 

C1008  Steel 

0.02 

1.2 

0.3 

100 

mm 

MAR  200-Steel 

0.01 

4.8 

1.2 

40 

-5.0 

Armco  Iron 

0.05 

2.0 

0.5 

70 

Tantalum 

0.003 

4.8 

1.2 

10 

-3.0 

1020  Steel 

0.01 

2.4 

0.6 

1 

■a 

MAR-250  Steel 

0.01 

3.2 

0.8 

50 

-3.0 

AF1410  Steel 

0.01 

5.6 

1.4 

20 

-4.0 

28 


Based  on  the  sensitivity  study  discussed  in  the  preceding  section,  the 
failure  model  parameters  are  estimated.  Out  of  the  six  parameters  (a^,  s, 
f  f^,  fi ,  and  N)  ,  determination  of  the  first  three  is  fairly  straightfor¬ 
ward.  The  calculated  value  for  a.,  was  16  kbar,  as  determined  from  the 
measured  value  of  AV^ .  In  order  to  minimize  the  number  of  model  parameters, 
by  taking  advantage  of  the  fact  that  s  and  f  are  less  sensitive  to  the 

failure  processes,  arbitrary  values  were  assigned  for  s  -  0.25  and  f  -  1. 

N  cr 

Also,  the  spall  signal  did  not  show  significant  differences  for  s  -  2,  4,  and 
8.  Similarly,  for  f  >  0.5,  the  results  showed  similar  trends;  therefore,  a 
theoretical  value  of  1  was  chosen  for  f 

cr 

The  remaining  parameters  /3,  N,  and  f^  were  determined  by  adjusting  them 
until  the  simulated  free  surface  velocity  matched  well  with  the  experimental 
data.  Following  the  guidelines  discussed  in  the  preceding  section,  f  was  ad¬ 
justed  to  approximately  match  the  spall  signal  arrival  time.  Then  f)  was 
modified  until  the  average  initial  slope  of  the  spall  signal  matched  the  ex¬ 
perimental  data.  Finally,  the  value  of  N  was  chosen  so  that  the  peak  velocity 
of  the  spall  signal  rebound  matched  with  the  data.  The  corresponding  failure 
model  parameters  for  OFHC  copper  are  given  below. 

a  s  f.  0  N  f 

N  1  cr 

(kbar)  (kbar) 

16  4  0.01  65  -2.4  1 

The  comparison  between  the  model  simulation  and  the  experimental 
velocity  history  data  for  OFHC  copper  is  shown  in  Figure  7.  Since  the  HEL  of 
OFHC  copper  is  extremely  low  compared  to  the  shock  stress,  the  BP  model 
predicted  HEL  (knee  at  A  in  Fig.  7)  did  not  match  well.  However,  for  copper, 
the  HEL  value  is  not  considered  critical  because  plastic  yielding  occurs  at 
(uniaxial)  stress  below  1  kbar.  Fine  tuning  of  the  BP  model  parameters  will 
certainly  improve  the  matching,  but  this  will  not  significantly  influence  t r. • 
failure  model  constants.  The  rounding  of  the  simulated  velocity  profile 
around  P  and  R  in  Figure  7  indicates  either  the  value  of  artificial  viscosi: 
coefficient  is  excessive  or  the  strain  rate  sensitivity  (n-0.4)  assumed  for 
copper  in  the  BP  model  is  large.  The  matching  beyond  the  point  S  is  control 
led  by  the  RDG  model  constants.  In  general,  the  overall  match  between  the 
simulation  and  the  data  seems  to  be  very  good.  This  demonstrates  the  abilit. 
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Fiqure  7.  Comparison  Between  Model  Simulation  and  the  Experimental  Data 
of  the  Free  Surface  Velocity  History. 


Co  calibrate  che  model  constants  using  plate  impact  test  data  as  well  as  the 
RDG  model  capability  to  accurately  reproduce  the  spall  signal. 

To  further  verify  the  model  constant  calibration  using  a  plate  impact 
test  data,  spall  in  HY100  steel  was  considered.  In  this  case,  a  manganin  gage 
recorded  stress  history  is  the  experimental  data.  Recall  that  the  failure 
model  constants  for  OFHC  copper  were  obtained  using  the  velocity  history.  The 
experimental  data  were  obtained  only  lor  a  time  duration  of  4  microseconds. 

In  this  experiment  the  release  waves  from  the  edges  of  the  flyer  plate  arrived 
at  the  gage  location  around  4.2  microseconds.  The  model  constants  were  ob¬ 
tained  by  matching  the  experimental  stress  signal  until  this  time.  The 
corresponding  match  is  shown  in  Figure  8a.  The  modeling  of  the  material  be¬ 
havior  before  any  void  nucleated  in  the  target  is  excellent.  This  indicates 
that  the  BP  model  constants  reproduced  the  HEL  data  well.  Also,  the  strain 
rate  sensitivity  of  HY100  steel  seemed  to  be  modeled  correctly  as  one  can  in¬ 
terpret  from  the  shape  of  the  simulated  signal  around  the  stress  peak. 
Rajendran  et  al.  [24]  modeled  the  SHB  test  accurately  using  the  same  BP  model 
constants  for  HY100  steel.  The  ability  to  model  both  the  SHB  and  plate  impact 
tests  makes  the  model  constants  relatively  more  general  than  by  matching  only 
one  of  these  two  test  configurations.  In  Figure  8b,  the  computed  stress  his¬ 
tory  at  the  spall  plane  is  shown.  The  compressive  stress  is  shown  as 
positive,  and  the  tensile  stress  is  shown  as  negative.  The  material  ex¬ 
periences  a  compressive  stress  of  around  lOGPa.  When  the  material  goes  into 
tension,  the  stress  peak  only  reaches  2.8  GPa.  This  reduction  in  stress  under 
tension  is  due  to  void  softening.  The  stress  relaxes  to  zero  as  the  material 
spalls.  Modeling  of  this  time  dependent  stress  relaxation  is  the  most  impor¬ 
tant  feature  of  the  RDG  model. 

The  failure  model  constants  were  accurately  determined  as  can  be  seen 
from  the  excellent  comparison  between  the  model  simulation  '■nd  the  experimen¬ 
tal  data  in  Figure  8.  Since  the  BP  model  constants  for  the  intact  material 
and  the  RDG  model  constants  for  the  void-containing  aggregate  have  been  ac¬ 
curately  determined,  the  entire  stress  history  could  be  reproduced  accurately 
Of  course,  this  statement  assumes  that  the  equations  of  state  for  the  flyer, 
target,  and  the  PMMA  are  accurately  modeled  in  the  computer  code  simulation. 
The  ability  of  the  model  to  reproduce  the  entire  spall  signal  clearly 
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demonstrates  that  the  nucleation  and  growth  models  follow  the  physical 
processes  realistically. 

RDG  model  constants  were  also  determined  for  1020,  C1008 ,  and  MAR- 200 
steels,  Armco  Iron,  and  pure  tantalum  by  matching  the  experimental  data  as 
shown  in  Figures  9  and  10.  The  ability  of  the  RDG  model  in  modeling  the  spall 
failure  is  truly  outstanding.  The  RDG  model  generated  spall  signal  for  tan¬ 
talum  matched  extremely  well  with  the  data  up  to  3.4  microseconds.  The 
difference  between  the  model  and  experiment  could  be  due  to  the  uncertainty  in 
the  experimental  data.  In  this  case  (#7-1302),  the  plate  impact  configuration 
was  such  that  the  spall  plane  location  was  very  close  to  the  manganin  gauge 
(about  one  millimeter) . 

In  the  simulations,  behavior  of  the  intact  material  can  also  be 
described  by  the  JC  model  instead  of  the  BP  model.  Using  the  JC  and  RDG 
models,  the  plate  impact  test  #1299  on  HY100  steel  was  simulated.  The  JC  con¬ 
stants  for  the  HY100  steel  were  obtained  by  Johnson  and  Holmquist  [34]  using 
torsional  SHB  test  data.  In  the  simulation,  the  authors  initially  used  the 
RDG  model  constants  for  HY100  that  were  determined  in  conjunction  with  the  BP 
model.  In  Figure  11  the  stress  histories  obtained  from  1)  test  data,  2) 
simulation  using  BP  model  and  the  RDG  model,  and  3)  simulation  using  JC  model 
and  the  RDG  model  are  shown.  The  matching  between  1  and  2  are  excellent  which 
is  not  too  surprising  because  the  RDG  model  constants  were  calibrated  in  con¬ 
junction  with  the  BP  model.  When  the  same  RDG  model  constants  were  employed 
in  the  simulation  in  conjunction  with  the  JC  model,  the  spall  signals  did  not 
match  exactly.  This  slight  difference  could  be  attributed  to  the  material 
characterizations  of  HY100  steel  using  different  experimental  configurations 
and  sources.  In  other  words,  the  descriptions  of  HY100  steel  by  the  BP  model 
and  the  JC  model  are  slightly  different;  therefore,  this  difference  will  in¬ 
troduce  slight  alterations  in  the  RDG  model  constants.  To  determine  the  RDG 
model  constants  that  correspond  to  the  use  of  JC  model  for  the  intact 
material,  plate  impact  simulations  were  made  and  the  constants  were  deter¬ 
mined.  A  slight  adjustment  to  the  f^  to  a  value  of  0.01  led  to  a  better  fit 
between  the  failure  model  and  the  test  data.  It  can  be  seen  from  Figure  12 
that  the  spall  signal  is  matched  extremely  well  by  both  BP  and  JC  models  with 
the  use  of  appropriate  RDG  model  constants. 
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An  important  aspect  of  the  modeling  is  the  model's  ability  to  reproduce 
an  experiment  which  was  not  used  for  determining  the  model  constants.  Also 
the  model's  ability  to  reproduce  a  failure  process  in  some  other  stress  -  strain 
configuration  is  important.  To  check  this  aspect,  a  new  experiment  (#1288)  on 
OFHC  copper  was  considered.  In  this  experiment,  the  stress  history  at  the  in¬ 
terface  of  the  target  and  the  back  plate  was  measured  using  a  manganin  gage. 
The  flyer  and  target  thicknesses  and  the  impact  velocity  were  all  different 
than  in  experiment  #538.  The  new  experiment  #1288  was  simulated  using  the 
same  model  constants  that  were  determined  using  experiment  #538  to  check  the 
model's  ability  to  predict  the  measured  stress  history.  In  Figure  13,  a  com¬ 
parison  between  the  model  prediction  and  the  data  is  shown.  The  model 
reproduced  the  measured  spall  signal  extremely  well  as  can  be  seen  from  this 
figure . 
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5.  MODELING  OF  VOID  COLLAPSE  UNDER  MULTIPLE  SHOCKING 


In  real  life  application  problems,  such  as  projectile  penetration  into 
a  target,  explosively  compacted  metal  operations,  etc,,  the  loading  conditions 
often  lead  to  multiple  shocking  in  the  material.  Under  these  conditions,  void 
growth  as  well  as  void  collapse  occur.  The  dynamic  properties  of  the  shocked 
material  are  modified  due  to  loading  history  effects.  To  further  evaluate  the 
RDG  model,  a  double  flyer  plate  impact  experiment  was  simulated  with  the  EPIC- 
2  finite  element  code  [29],  Previously  calibrated  RDG  model  parameters  were 
used  in  the  simulation.  The  computed  stress  history  was  compared  with  the  ex¬ 
perimentally  measured  stress  history.  A  simple  void  collapse  criterion,  based 
on  a  critical  void  volume  fraction,  was  employed  to  obtain  a  better  fit  to  the 
experimental  data. 

5.1  Double  Flyer  Plate  Impact  Experiment 

Yaziv  [36]  and  Yaziv  et  al .  [37]  introduced  a  double  flyer  plate  impact 
technique  for  examining  the  dynamic  properties  of  shock  damaged  materials . 

This  technique  differs  from  a  conventional  plate  impact  experiment  in  that  two 
flyer  plates  are  used,  separated  by  a  small  gap,  as  shown  schematically  in 
Figure  14. 

The  first  flyer  has  a  lower  shock  impedance  than  the  second  flyer. 

When  the  first  flyer  impacts  the  target  plate,  microvoid  nucleation  and  growth 
occur  in  a  localized  area  of  maximum  tensile  stress  in  the  target.  The 
reshocking  by  the  second  flyer  reverses  the  damage  by  closing  the  microvoids. 

Grove  et  al.  [38]  simulated  an  experiment  (#829)  in  which  the  first 
flyer  was  a  2  mm  thick  aluminum  plate  and  the  second  flyer  was  a  3  mm  thick 
copper  plate.  The  spacing  between  the  flyers  was  0.25  mm,  and  the  impact 
velocity  was  approximately  336  m/s.  The  target  was  a  4  ram  thick  OFHC  copper 
plate  with  a  FMMA  backing.  The  stress  history  at  the  interface  of  the  target 
and  PMMA  was  measured  with  a  manganin  stress  gauge. 

The  double  flyer  plate  impact  experiment  was  simulated  using  the  EPIC 
finite  element  code  [29],  EPIC-2  material  models  were  used  to  describe  both 
flyers  and  the  PMMA.  The  PMMA  was  modeled  as  a  rate  independent  elastic - 
perfectly  plastic  material  with  a  yield  stress  of  0.18  GPa.  The  RDG  model  w.. 
used  to  describe  the  plastic  flow  and  dilatation  in  the  OFHC  copper  target 
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Figure  12.  Comparisons  between  the  spall  signals  generated  using 
JC  and  BP  models  for  different  values  of  RDG  model 
constants . 


oe 


(Bd3)  ssaj^s 


40 


Figure  13.  Comparison  between  model  prediction  and  the  measured 
stress  history  (OFHC  copper) . 


plate.  The  previously  determined  RDG  model  parameters  for  OFHC  copper  had 
been  used  in  this  simulation. 

In  the  first  simulation  of  the  double  flyer  plate  experiment,  the  RDG 
model  predicted  the  initial  spall  signal  reasonably  well  but  was  unable  to  ac 
curately  reproduce  the  experimental  recompaction  signal,  as  indicated  by 
Figure  15.  The  plastic  wave  slope  of  the  simulated  reshocking  exhibited  ramp 
ing  with  a  long  rise  time  indicating  that  the  predicted  pore  collapse, 
described  by  the  same  equations  as  void  growth,  was  unrealistically  slow  and 
required  modification.  Grove  et  al.  [38]  introduced  a  complete  collapse  of 
the  microvoids  in  a  particular  element  if  that  element's  dilatation  rate  was 
negative  and  its  void  volume  fraction  was  less  than  or  equal  to  some  critical 
value.  Figure  16  shows  the  excellent  fit  to  the  experiment  obtained  when  the 
critical  void  volume  fraction  for  complete  collapse  was  five  percent. 

The  double  flyer  plate  experiment  without  void  growth  was  also  simu¬ 
lated  so  that  the  microvoids  effect  on  the  reshock  signal  could  be 
investigated.  Figure  17  compares  the  two  simulations,  with  and  without  void 
growth.  It  is  interesting  to  note  that  the  dip  at  point  F  occurs  at  the  same 
time  in  both  curves,  indicating  that  this  feature  is  directly  related  to  the 
experimental  configuration.  Figure  18  shows  an  x-t  (distance-time)  diagram 
for  the  double  flyer  plate  impact  experiment,  assuming  no  void  growth.  In 
this  figure,  the  solid  and  dashed  lines  represent  shock  and  release  waves, 
respectively.  The  release  and  shock  waves  that  arrive  at  points  E  and  F  are 
solely  due  to  the  impedance  mismatch  between  the  first  flyer  and  the  target. 
Because  the  first  flyer  has  a  lower  impedance  than  the  target,  it  separates 
from  the  target  at  about  0.75  microseconds.  As  Figure  18  indicates,  the  in¬ 
itial  shock  wave  from  the  second  flyer  reflects  off  the  free  surface  of  the 
first  flyer  as  a  release  wave  because  of  the  gap  formed  between  the  first 
flyer  and  the  target.  This  release  wave  can  be  traced  to  point  E.  When  the 
first  flyer  impacts  the  target  again,  another  shock  wave  is  produced,  and  thi 
shock  wave  can  be  traced  to  point  F  in  the  x-t  diagram. 

In  the  absence  of  microvoid  nucleation  and  growth,  the  shock  wave  from 
the  second  flyer  would  arrive  at  point  C.  When  void  growth  occurs,  however, 
this  signal  is  delayed  because  of  the  reduced  wave  speed  in  the  porous  region 
of  the  target.  Arrival  of  this  signal  at  point  C'  (in  Figure  17)  signifies 
that  the  microvoids  have  completely  collapsed.  Since  the  release  and  shock 
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Figure  14.  Experiment  Configuration  for  Double 
Flyer  Impact  Experiment. 
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Figure  15.  Stress  history  comparison  for  originial  RDG  model. 
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Figure  16.  Stress  history  comparison  for  modified  RDG  model, 
with  a  critical  void  volume  fraction  of  5%  for 
complete  collapse. 


waves  ac  points  E  and  F  arrive  after  point  C' ,  it  appears  that  they  are  unaf¬ 
fected  by  the  void  growth  phenomenon. 

In  summary,  a  double  flyer  plate  impact  experiment  was  successfully 
simulated  using  the  RDG  model.  Using  previously  calibrated  constants,  the  RDG 
model  was  able  to  reproduce  the  initial  spall  signal,  but  the  reshock  signal 
did  not  compare  as  well  with  the  experimental  stress  history.  The  RDG  model 
was  modified  to  include  a  simple  void  collapse  criterion  in  which  the 
microvoids  were  forced  to  completely  collapse  below  a  critical  void  volume 
fraction  when  the  dilatation  rate  was  negative.  Microvoid  collapse  in  OFHC 
copper  was  successfully  modeled  using  precalibrated  RDG  model  constants  and  a 
critical  void  volume  fraction  for  complete  collapse  of  5  percent.  Thus,  the 
RDG  model  has  been  extended  to  describe  microvoid  collapse  as  well  as  a  his¬ 
tory  dependent  failure  process. 
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Figure  17.  Comparison  of  simulated  stress  histories  from  double 
flyer  impact  configuration,  with  and  without  void 
growth . 
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Figure  18.  x-t  diagram  for  the  double  flyer  impact  experiment, 
assuming  no  void  growth.  This  corresponds  to  the 
dashed  line  in  Figure  17. 
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6.  RDG  MODEL  APPLICATIONS 


The  dynamic  ductile  failure  also  occurs  in  other  types  of  configura¬ 
tions.  Tensile  necking  in  a  split  Hopkinson  bar  test  and  spallation  in  a 
target  plate  impacted  by  a  long  rod  penetrator  are  a  few  examples.  These  two 
configurations  are  investigated  in  the  following  sections. 

6.1  Tensile  Necking 

To  further  demonstrate  the  capability  of  RDG  model,  several  two  dimen¬ 
sional  axi- symmetric  configurations  were  considered.  The  first  case 
considered  the  problem  of  dynamic  tensile  necking.  The  main  objective  of  this 
two  dimensional  analysis  was  to  examine  the  stress  states,  evolution  of  neck¬ 
ing,  void  growth  as  influenced  by  the  plastic  deformation,  and  ductility.  To 
initiate  necking,  a  shallow-notched  copper  specimen  was  considered.  This  is 
shown  schematically  in  Figure  19.  The  length  and  uniform  diameter  of  the 
specimen  are  0.35  and  0.125  inches  respectively.  The  minimum  diameter  at  the 
notch  is  0.105  inches.  The  length  represents  the  actual  gauge  length  of  a 
standard  SHB  tensile  specimen. 

The  shallow  notched  specimen  was  modeled  using  the  quadrilateral  ele¬ 
ment  option  in  the  EPIC-2  code.  Due  to  specimen  symmetry,  it  is  sufficient  to 
model  only  one  quarter  of  the  specimen  as  shown  in  Figure  20.  The  Bodner- 
Partom  (BP)  model  [22]  was  used  to  describe  the  high  strain  rate  behavior  of 
the  copper  specimen.  The  BP  model  constants  for  copper  are  given  in  Table  1 
A  typical  experimentally  measured  SHB  pull  velocity  was  applied  to  the 
specimen  in  the  simulation.  The  velocity  was  linearly  increased  from  zero  to 
50.8  m/sec  (2000  inches/sec)  for  40  microseconds,  and  afterwards,  it  was  kept 
constant . 

The  RDG  model  [1]  was  used  to  describe  the  necking  evolution,  and  the 
corresponding  constants  are  given  in  Table  4.  A  mean-effective  plastic  strai: 
based  Gaussian  distribution  as  described  by  equation  [9]  was  used  for  void 
nucleation.  The  corresponding  RDG  model  constants  are  given  in  Table  4. 

Since  the  stress  level  under  uniaxial  stress  state  is  considerably  lower  tha:. 
the  level  under  uniaxial  strain  condition  (plate  impact  test),  the  void 
nucleation  was  assumed  to  occur  entirely  due  to  large  plastic  flow  around  thu 
inclusions  and  oxide  particles.  In  a  uniaxial  tensile  test,  the  necking  ofti 
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Figure  20.  Finite  element  mesh  for  shallow-notched  SHB  specimen. 

The  velocity  loading  history  is  shown  in  the  inset. 
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progresses  under  significant  plastic  flow  and  a  moderate  triaxial  stress 
state. 

In  the  analysis,  the  time  histories  of  effective  stress,  effective 
plastic  strain,  void  volume  fraction,  and  triaxiality  of  the  stress  state 
(ratio  of  pressure  and  effective  stress)  were  evaluated  for  elements  #82, 

#182,  #282,  and  #382  as  shown  in  Figure  20.  Variation  of  pressure,  stress 
components,  strength,  and  other  variables  along  a  radius  (snapshots)  for  dif¬ 
ferent  time  intervals  have  also  been  considered  in  the  analysis. 

The  stress  triaxiality  parameter  (P/Y)  contours  are  plotted  in  Figure 
21.  While  this  parameter  is  increasing  at  the  specimen  center,  the  stress 
state  is  under  uniaxial  stress  condition  (P/Y  is  around  1/3)  for  the  uniform 
shaded  region.  The  triaxiality  ratio  becomes  negative  (compressive)  around 
region  A  as  can  be  seen  from  Figure  21b.  The  unloading  of  tne  uniform  section 
(region  B)  and  the  compressive  loading  at  A  make  the  interpretation  of  the 
triaxiality  contours  complicated.  In  the  presence  of  increased  triaxiality, 
the  void  growth  also  increases  as  shown  in  Figure  22.  The  void  volume  frac¬ 
tion  contours  at  t-0,  20,  50,  and  60  microseconds  are  shown.  Initially,  a 
small  amount  of  voids  nucleate  near  the  stress -free  notch  surface.  The  void 
distribution  spreads  toward  the  specimen  center  with  increasing  necking.  By 
about  50  microseconds,  the  voic.  ontent  at  the  center  region  increases 
rapidly  to  3  to  4  percent.  Due  to  increasing  triaxial  tensile  state,  rapid 
void  growth  occurs  during  the  next  10  microseconds;  by  t-60 ,  the  necked  region 
is  filled  with  voids,  over  40  percent  located  at  the  center  (see  Figure  22d) . 
Although  the  void  nucleated  near  the  surface,  the  failure  initiation  even¬ 
tually  occurs  at  the  center  of  the  specimen.  Hancock  and  Brown  [39]  provided 
the  physical  evidence  to  support  failure  initiation  at  the  center  in  notched 
tensile  specimens  for  various  notch  geometries  under  quasi-static  loading  con¬ 
ditions  . 

The  time  history  plots  of  effective  plastic  strain  and  void  volume 
fraction  for  regions  near  the  neck  (local,  element  #382)  and  away  from  the 
neck  (uniform,  element  #82)  are  given  in  Figure  23.  Initially,  the  plastic 
strains  are  identical  at  both  regions.  However,  the  local  strain  at  the  cen¬ 
ter  of  the  specimen  starts  deviating  from  the  uniform  strain  slowly  due  to 
enhanced  void  growth.  The  void  nucleation  process  occurs  during  the  first  40 
microseconds,  and  rapid  growth  initiates  around  t-50.  Beyond  this  time,  while 
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Figure  23.  Time  histories  of  effective  plastic  strain  and  void 

volume  fraction  for  elements  #82  (uniform  section)  and 
#382  (local  section. 


the  local  strain  continues  increase  due  to  void  growth,  the  uniform  strain  at 
element  #82  (see  Figure  20)  reaches  a  maximum  value  of  20  percent  and  stopped 
increasing.  This  is  obviously  due  to  unloading  of  the  sections  away  from  the 
local  necking  region,  shown  in  Figures  24-26. 

The  stress -time  history  for  the  local  and  uniform  regions  are  shown  in 
Figure  24.  Also  shown  in  this  figure  is  the  void  volume  fraction  history  for 
elements  #82  and  #382.  Since  these  two  elements  are  very  close  to  the  axis  of 
symmetry,  the  hoop  and  radial  stresses  are  identical  (Figure  24).  For  the 
uniform  element  (#82) ,  the  hoop  and  radial  stresses  are  zero  as  one  would  ex¬ 
pect.  However,  the  stress  state  in  the  local  element  #382  begins  almost 
uniaxial  and  becomes  triaxial  due  to  necking.  The  hoop  and  radial  stresses 
increase  as  the  void  growth  occurs  in  the  local  region  generating  a  high  tri¬ 
axial  stress  state  locally.  Between  40  and  50  microseconds,  this  stress  state 
accelerates  the  necking  process.  As  rapid  void  growth  occurs  at  50 
microseconds,  the  axial  as  well  as  hoop  and  radial  stresses  relax  (at  point 
A')  in  local  element  #382.  The  word  'relax'  indicates  the  decreasing  stresses 
with  increasing  plastic  strains.  However,  the  axial  stress  actually  unloads 
in  uniform  element  #82  (at  point  A  in  Figure  24)  due  to  rapid  localization  of 
the  deformation.  Since  no  void  growth  occurs  in  this  uniform  region,  the 
decreasing  stresses  are  only  due  to  unloading.  In  regions  away  from  the  neck, 
the  material  is  intact  without  any  voids,  and,  therefore,  the  possibility  of 
any  strength  degradation  leading  to  stress  relaxation  does  not  exist.  To  fur¬ 
ther  demonstrate  these  aspects,  effective  stress  versus  effective  plastic 
strain  is  shown  in  Figure  25.  The  void  volume  fraction  history  has  also  been 
included.  The  stress  (strength  degradation)  relaxation  beyond  point  A  can  be 
seen  in  this  figure.  The  plastic  strain  continually  increases  because  of  the 
void  growth  induced  plastic  flow.  As  the  void  volume  fraction  increases,  the 
strength  continues  to  decrease. 

The  mean  stress  (pressure)  histories  at  various  locations  are  plotted 
in  Figure  26.  As  mentioned  earlier,  #82  represents  the  uniform  region  and 
#382  represents  the  local  region.  The  elements  #282  and  #182  are  the  inter¬ 
mediate  regions  (see  Figure  20) .  The  stress  triaxiality  (defined  as  the  ratio 
of  mean  stress  and  effective  stress)  due  to  the  initial  notch  leads  to  a 
higher  mean  stress  level  in  the  local  region  before  the  occurrence  of  necking. 
When  necking  initiates  around  40  microseconds,  the  regions  (#82  and  #182)  away 
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Figure  24.  Time  history  plots  for  the  stress  components  and  void 
volume  fraction. 
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Figure  27.  Measure  of  triaxiality  with  respect  to  time  for  the 
uniform  (#82)  and  local  (#382)  sections. 

The  dashed  line  represents  the  uniaxial  stress  state 


from  the  neck  start  unloading.  While  the  pressure  decreases  because  of  un¬ 
loading  in  these  regions,  the  mean  stress  relaxes  due  to  void  growth  in  the 
local  regions. 

The  stress  triaxiality  is  plotted  with  respect  to  time  for  elements  #82 
and  #382  in  Figure  27.  The  uniform  region  remained  closely  under  uniaxial 
stress  state.  The  ratio  of  mean  stress  and  effective  stress  was  almost  1/3, 
as  one  would  expect.  The  slight  deviation  from  this  value  before  20 
microseconds  is  numerical,  and  later  this  ratio  remains  close  to  1/3.  For  the 
element  in  the  local  region,  the  stress  triaxiality  is  initially  greater  than 
1/3  due  to  the  notch.  However,  as  necking  occurs  between  40  and  50 
microseconds,  the  local  triaxiality  rapidly  increases  in  Figure  27.  Beyond 
point  B,  rapid  void  growth  leads  to  strength  degradation  and  failure  initia¬ 
tion.  For  further  understanding  the  loading  conditions,  the  loading  paths  are 
shown  at  different  locations  with  respect  to  the  pressure  dependent  yield  sur¬ 
faces  in  Figure  28. 

The  loading  paths  for  elements  #82,  #182,  #282,  and  #382  are  described 

by  plotting  the  ratio  of  effective  stresses  of  the  aggregate  and  the  fully 

dense  matrix  (a  __/Y  )  and  the  ratio  of  mean  stress  and  effective  stress  of 
eff  m 

the  matrix  (P/Y^) .  Since  the  material  is  initially  void-free,  the  plastic 
yielding  is  pressure  independent.  Therefore,  the  loading  path  starts  at 
cr  ^/Y  -  1  as  shown  in  Figure  28.  The  dotted  lines  represent  the  elliptical 

yield  surfaces  for  various  levels  of  porosity.  In  element  #82,  a  small  amount 
of  void  nucleates  initially  and  grows  very  little.  The  corresponding  loading 
path  is  shown  between  points  A  and  B.  The  loading  remains  on  the  yield  sur¬ 
face  for  void  volume  fractions  less  than  0.3  percent.  In  a  similar  manner  the 
loading  path  for  element  #182,  which  is  away  from  the  notch  zones,  is  confined 
to  a  yield  surface  corresponding  to  f  -  0.01.  Because  of  enhanced  void  growth 
in  the  local  regions,  a  significant  amount  of  void  growth  occurred,  especially 
in  element  #382  as  between  points  C  and  D  in  Figure  28.  As  the  loading  path 
proceeds  along  the  shrinking  yield  surfaces,  the  material  strength  degrades  to 
less  than  80  to  90  percent  of  the  original  intact  material,  indicating  failure 
of  the  material. 

For  completing  the  discussions  of  the  RDG  model  capabilities,  addi¬ 
tional  analyses  in  terms  of  snap  shots  and  additional  contour  plots  are 
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presented.  The  snap  shots  are  plots  that  describe  the  distribution  of  a  vari¬ 
able  with  respect  to  position  (for  example,  radial  distance  from  the  axis-of- 
symmetry) .  In  the  snap  shot  plots,  the  radius  is  normalized  with  respect  to 
the  initial  notch  radius.  In  Figure  29,  the  snap  shots  of  void  volume  frac¬ 
tion  for  the  three  different  z  positions,  shown  earlier  in  Figure  20,  are 
plotted.  The  z  position  describes  the  position  along  the  specimen  length  that 
is  away  from  the  minimum  cross  section  (c/s)  of  the  shallow  notched  tensile 
specimen.  Therefore,  z-0  represents  the  minimum  cross  section.  Since  the 
triaxiality  is  largest  at  the  center  of  the  specimen  (r-0  and  z-0) ,  the  void 
volume  fraction  is  maximum  at  the  center  and  decreases  toward  the  free  surface 
(r/r  -1).  The  present  analysis  is  for  a  dynamically  stretched  tensile 
specimen.  The  loading  duration  is  of  the  order  of  microseconds. 

Interestingly,  the  numerical  results  are  very  similar  to  a  statically  deformed 
specimen.  These  analyses  compared  qualitatively  well  with  the  results 
reported  by  Norris  et  al.  [4-0].  They  analyzed  quasi-static  necking  deforma¬ 
tion  using  a  dynamic  finite  difference  computer  program.  However,  in  the 
present  work,  the  problem  is  considered  as  inherently  dynamic  and  has  been 
treated  as  a  shock  wave  propagation  problem.  Due  to  several  reverberations  of 
the  shock  waves,  th*3  deformation  is  homogeneous,  and  the  problem  is  dominated 
by  inertial  loading  rather  than  the  shock  wave  effects. 

The  snap  shot  of  normalized  effective  stress  of  the  voided-aggregate  is 
given  in  Figure  30.  At  t-20  microseconds,  the  effective  stresses  of  the 
matrix  and  the  aggregate  are  almost  the  same  due  to  insignificant  void  growth; 
therefore,  the  ratio  is  1.  When  the  void  growth  occurs  at  z-0  c/s,  the  ag¬ 
gregate  strength  degrades  and  the  ratio  becomes  less  than  one.  At  40 
microseconds,  the  c/s  at  z-0 . 84  mm  and  the  minimum  c/s  develop  significant 
voids.  The  z-0  c/s  degrades  faster  than  the  section  at  0.84  mm.  Since  the 
void  volume  fraction  is  initially  evenly  distributed  across  the  c/s  (not  shown 
in  the  figure),  the  reduction  in  strength  is  almost  uniform.  However,  later 
at  50  microseconds,  the  non-uniform  void  growth  in  the  intermediate  section 
leads  to  rapid  strength  degradation  near  the  axis  (r/ag-0)  as  shown  between  A 
and  B  in  Figure  30.  However,  because  of  a  lack  of  stress  triaxiality  near  the 
surface,  the  void  growth  was  not  significant;  correspondingly,  the  strength 
did  not  significantly  degrade.  At  the  minimum  c/s,  due  to  the  enhanced  void 
growth  and  void  coalescence  across  most  of  the  c/s,  the  strength  degraded  nor 
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Figure  29.  Radial  distribution  of  void  volume  fraction  at  these 
different  positions  along  the  specimen  length. 


and 


Figure  30.  Normalized  effective  stress  la  y  )  distribution 

etf'  ■‘m 

along  the  radius  for  different  z  positions  and  time 
intervals . 


than  90  percent.  This  indicates  failure  initiation  leading  to  total  separa¬ 
tion  of  the  material. 

For  comparison,  the  necking  process  was  also  simulated  without  the 
voids.  This  simulation  will  be  referred  to  as  case  2  in  the  following  dis¬ 
cussions.  The  Bodner-Partom  model  described  the  material  behavior.  The 
solution  was  carried  out  for  60  microseconds.  Since  the  boundary  conditions 
for  the  two  simulations  were  the  same,  the  top  section  was  moved  (pulled) 
identically.  Therefore,  at  any  given  instant,  the  specimen  length  was  the 
same  for  both  cases.  Plastic  strains  were  plotted  based  on:  the  ratio  of 
current  length  and  original  length;  uniform  strain  in  element  #82  for  void- 
free  material  (case  2);  and,  uniform  strain  in  element  #82  for  void-containing 
material  (case  1)  in  Figure  31.  As  one  would  expect,  the  average  strain  based 
on  the  overall  specimen  length  was  larger  than  the  uniform  strains  in  element 
#82  for  the  two  cases.  While  in  case  2  the  uniform  section  continues  to 
deform  at  60  microseconds,  deformation  no  longer  occurs  in  case  2  due  to  void 
growth  controlled  necking  process.  For  this  reason,  the  ductility  reduced 
significantly  in  the  presence  of  void  growth. 

The  effective  stress  versus  effective  plastic  strain  in  the  local  ele¬ 
ment  for  the  cases  with  and  without  voids  is  shown  in  Figure  32.  In  the  void- 
free  case,  the  necking  process  continues  beyond  60  microseconds.  Since 
m.. serial  degradation  does  not  occur  in  this  case,  the  effective  stress 
(strength)  follows  the  true  stress-  true  strain  curves  generated  by  the  BP 
model  for  OFHC  copper  (Figure  32).  The  drooping  portion  of  the  stress-strain 
curve  (between  points  A  and  B)  for  the  void- growth  influenced  necking  case  is 
due  to  strength  degradation.  Even  though  the  loading  (increasing  plastic 
strain)  increases  beyond  point  A,  the  stress  continually  drops  indicating 
void- softening  in  the  material. 

The  effective  plastic  strain  variations  along  a  radius  for  the  necking 
process  with  and  without  voids  are  shown  in  Figure  33.  The  variation  along 
the  minimum  radius  (local)  indicates  that  the  strain  reached  a  maximum  near 
the  surface  and  not  at  the  axis  for  both  cases.  These  results  differ  from  the 
tensile  simulation  results  reported  by  Norris  et  al.  [40]  in  which  the  plastic 
strain  reached  maximum  at  the  axis.  They  considered  a  quasi-static  necking 
problem  in  which  the  inertial  effects  were  neglected.  Since  the  present  neck¬ 
ing  analysis  considers  the  inertial  effects,  the  difference  in  the  effective 
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Figure  31.  Comparisons  of  failure  strain  for  tensile  necking 
with  and  without  voids.  The  solid  line  without 
symbols  is  an  average  measure  of  strain  based  on  the 
pull  veloc  Ity . 
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plastic  strain  distribution  is  attributed  to  the  inertial  loading.  Also,  the 
stress-strain  response  was  modeled  as  an  isothermal  process.  In  the  simula¬ 
tions,  the  effects  due  to  the  adiabatic  process  are  not  considered. 

The  local  plastic  strain  levels  in  the  void-growth  case  are  larger  than 
that  in  the  void- free  case.  The  obvious  reason  is  that  void  growth  and 
pressure -dependent  yield  enhance  the  plastic  flow.  While  loading  continues  in 
the  local  region,  the  uniform  section  unloads  because  of  enhanced  necking 
process  due  to  void  growth.  Since  the  necking  was  not  very  severe  in  the 
void- free  case,  the  loading  continues  in  the  uniform  section  as  shown  earlier 
in  Figure  32.  The  effective  strain  along  the  uniform  cross  section  continued 
to  increase  in  this  case. 

The  snapshots  of  stress  triaxiality  for  the  two  cases  were  also 
plotted.  The  triaxiality  was  maximum  at  the  local  (minimum)  section  with  void 
growth,  especially  near  the  axis  (r-0)  as  can  be  seen  from  Figure  34.  Other 
investigators  [40-41]  have  observed  similar  behavior  under  quasi-static  load¬ 
ing  conditions  with  and  without  void  growth.  The  triaxiality  levels  are 
relatively  low  near  the  surface,  as  one  would  expect.  In  the  void- free  case, 
the  stress  triaxiality  is  equal  to  0.3  in  the  uniform  section.  However,  for 
the  uniform  section  of  the  void-growth  case,  the  triaxiality  near  the  axis  is 
actually  lower  than  0.3.  The  reason  for  this  can  be  explained  from  Figure 
35. 

Figure  35  shows  the  mean  stress  contour  plots  for  the  void-growth  case. 
The  material  under  compression  is  shown  by  the  shaded  region.  The  contour  was 
plotted  for  the  final  solution  time  (60  microseconds).  Failure  initiation  has 
occurred  due  to  large  void  content  near  the  central  portion  of  the  tensile 
specimen,  failure  initiation  has  occurred.  This  resulted  in  unloading  of  the 
uniform  section  and  further  loading  under  compression.  In  the  void- free  case 
(not  shown  in  the  figure) ,  the  uniform  section  unloaded  completely  but  did  not 
go  into  compression.  If  calculations  were  carried  out  beyond  60  microseconds, 
compressive  loading  might  also  occur  in  the  void- free  case.  However,  Norris 
et  al.  [40]  in  their  numerical  analysis  of  quasi-static  tensile  necking  ob¬ 
served  compressive  zones  away  from  the  local  region. 

The  deformed  configurations  of  the  shallow  notched  tensile  specimen  for 
the  cases  with  and  without  void  growth  are  shown  in  Figure  36.  Only  one  half 
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Figure  35.  Pressure  contour  plot  for  void-growth  influenced  necking  process 
The  shaded  region  is  under  tensile  loading,  and  the  rest  of  the 
region  experiences  compressive  loading. 


of  the  specimen  is  shown  for  each  case  for  comparison.  Also  included  was  the 
initial  geometry,  shown  by  the  dotted  line.  In  the  void-free  case,  the  neck 
profile  retained  the  original  notch  shape  while  in  the  void- growth  case  the 
profile  did  not  follow  the  notch  shape  as  can  be  seen  from  Figure  36. 

In  the  present  application,  the  material  was  modeled  as  strain  rate  de¬ 
pendent  and  highly  strain  hardening.  Additional  exercises  with  different 
material  behaviors  such  as  perfectly  plastic,  strain  rate  independent,  and 
thermal  softening,  with  and  without  voius  will  shed  additional  light  on  the 
evolution  of  necking. 

6.2  Application  of  RDG  Model  to  Rod  Penetration  Problem 

The  main  objective  of  this  simulation  is  to  show  the  spall  evolution  in 
the  target  using  the  RDG  model.  Until  now,  the  model  has  been  successfully 
used  to  describe  the  spallation  phenomenon  in  plate  impact  configuration  and 
void-growth  enhanced  tensile  necking  problems. 

An  EPIC-2  code  simulation  of  a  2-inch  long  and  1-inch  diameter  solid 
copper  rod  impacting  a  6 -inch  diameter  and  1  inch  thick  HY100  steel  at  an  im¬ 
pact  velocity  of  8333  ft/sec  was  considered.  The  corresponding  finite  element 
mesh  is  shown  in  Figure  37.  The  behavior  of  intact  copper  was  modeled  using 
the  Johnson-Cook  (JC)  model  [23]  and  the  HY100  steel  through  the  Bodner-Partom 
(BP)  model  [22].  The  constants  for  both  the  BP  and  JC  and  the  RDG  models  are 
given  in  Tables  1,  2,  and  4  respectively. 

The  numerical  solution  was  obtained  for  20  microseconds.  For  analysis, 
elements  near  the  axis  of  5yraraetry  of  the  target  and  at  a  distance  of  0.8125 
inches  from  the  impact  point  were  examined.  The  time  history  of  pressure, 
stress  components,  effective  stress,  and  damage  (void  volume  fraction)  were 
stored  for  post-processing.  Also  plotted  was  the  variation  of  these  variable.' 
along  the  target  radius  at  different  depths  for  various  time  intervals. 

Damage  contours  at  different  time  intervals  were  plotted  to  show  spall  cones 
in  the  target. 

In  Figure  38,  the  normalized  pressure  (mean  stress  /  aggregate 
strength)  is  plotted  with  respect  to  time.  The  evolution  of  damage  in  terms 
of  void  volume  fraction  is  also  shown.  The  normalized  pressure  it  a  measure 
of  triaxiality  of  the  'tress  state.  When  the  ratio  is  negative,  it  indicates 
that  the  material  element  is  under  compression,  and  void  growth  cannot  occur. 
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Figure  38.  Measure  of  Traixiality,  Ratio  of  Pressure  (mean  stress) 

and  Strength  (effective  stress),  and  Void  Volume  Fraction 
with  Respect  to  Time. 


To  generate  voids  in  an  intact  ductile  metal,  high  triaxial  stress  state  under 
tension  is  required.  The  derivations  cf  Rice  and  Tracey  [7]  and  McClintock 
[6],  and  the  quasi-static  experiments  of  Hancock  and  Mackenzie  [9]  supported 
the  concept  that  high  triaxial  tensile  state  enhances  void  nucleation  and 
growth.  The  RDG  model  which  is  also  based  on  this  concept  was  able  to  model 
the  void  growth  accordingly  as  can  be  seen  from  Figure  38.  In  a  tensile  test, 
before  the  occurrence  of  necking,  the  triaxiality  ratio  is  1/3.  However,  in  a 
penetration  experiment  configuration,  the  triaxiality  could  be  large  enough  to 
generate  spallation  in  the  target.  In  this  particular  simulation,  the  tensile 
triaxiality  reached  extremely  high  (>10  in  Figure  38). 

In  element  #736,  the  release  wave  from  the  edge  (point  A  in  Figure  37) 
arrives  at  about  4  microseconds  and  releases  the  compressive  pressure.  Note 
that  the  distance  of  this  element  from  the  edge  A  in  Figure  37  is  around  0.9 
inches.  The  wave  speed  in  HY100  steel  is  0.23  inches/microsecond. 

The  arrival  time  of  the  release  wave  from  the  stress  free  back  surface 
of  the  target  is  about  5  microseconds.  Note  that  the  initial  shock  wave  will 
reflect  as  a  tensile  wave.  This  tensile  wave  puts  the  element  under  a  high 
triaxial  stress  state  leading  to  void  nucleation  and  growth.  The  void  volume 
fraction  increases  while  the  tensile  stress  state  prevails.  As  the  material 
degrades  due  to  voids,  the  pressure  as  well  as  the  strength  relaxes  as  shown 
in  Figure  39. 

The  matrix  material  (intact  HY100  steel)  and  porous  aggregate  (HY100 
steel  with  porosity)  strengths  are  also  shown  in  this  figure  with  respect  to 
time.  For  further  understanding  and  evaluation  of  the  RDG  model  capabilities, 
the  evolution  of  void  volume  fraction  (damage)  has  also  been  included  in  the 
plot.  Initially,  the  material  (target  element)  is  intact;  therefore,  the  ag 
gregate  and  matrix  strengths  are  the  same,  and  the  curves  between  points  A  and 
ft  are  identical.  Void  nucleation  occurs  around  t  -  4.9  microseconds.  The 
void  generation  degrades  the  aggregate  strength  as  shown  by  the  sudden  drop  ot 
the  effective  stress,  Y^  to  the  point  X  and  remained  lower  compared  to  the 
matrix  strength  (as  shown  by  the  dotted  line).  The  matrix  material  strength 
(the  strength  of  the  undamaged  HY100)  obviously  remained  larger  than  the 
degraded  material  strength  as  one  can  see  from  Figure  39. 
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In  Figure  40,  the  loading  path  in  element  #736  was  plotted  with  respect 
to  the  pressure  dependent  yield  surfaces.  Yield  surfaces  for  various  values 
of  the  void  content  are  shown  in  this  figure.  These  surfaces  are  generated 
from  the  relationship  (Equation  16)  between  the  normalized  pressure  and  nor¬ 
malized  effective  strength.  The  pressure  and  the  effective  stress  are 
normalized  with  respect  to  the  matrix  flow  strength.  Initially,  the  material 
is  void  free,  and  therefore  the  effective  stress  (strength)  of  the  aggregate 
and  the  matrix  are  the  same.  For  this  reason,  the  loading  path  starts  at  the 
point  A  where  the  value  of  the  normalized  effective  stress  is  one.  When  small 
amounts  (<0.005)  of  voids  nucleate  between  point  A  and  B,  the  pressure  depend¬ 
ence  of  the  yielding  makes  the  aggregate  flow  strength  a  function  of  pressure. 
The  void-containing  aggregate  strength  is  lower  than  the  matrix  strength,  and 
the  ratio  between  these  two  strengths  becomes  less  than  one.  As  the  loading 
progresses,  rapid  void  growth  occurs  between  the  points  B  and  D  in  Figure  40. 
The  void  content  in  element  #736  increases  to  about  0.02  at  C.  The  flow 
strength  is  about  60  percent  of  the  original  intact  material  (matrix) 
strength.  Further  void  growth  between  points  C  and  D  degrades  the  flow 
strength  to  about  10  percent  of  its  original  value.  Beyond  point  D,  the 
material  completely  loses  its  load  carrying  capacity  in  tension,  and  the  flow 
strength  as  well  as  the  mean  stress  (pressure)  relaxes  to  zero. 

The  RDG  model  does  not  include  any  void  coalescence  model;  however,  the 
void  growth  model  has  been  formulated  in  a  manner  that  rapid  void  growth  oc¬ 
curs  at  a  certain  void  volume  fraction  level  indicating  void  coalescence. 

This  can  be  clearly  seen  in  Figure  40  between  points  D  and  F. 

The  extent  of  voids  generated  in  the  target  due  to  high  triaxial  ten¬ 
sile  stresses  are  shown  in  Figure  41.  The  void  distribution  along  the  radius 
at  element  #736  for  two  different  time  intervals  (5  and  20  microseconds ) are 
shown.  A  small  amount  of  voids  has  nucleated  and  grown  at  the  center  of  the 
specimen  between  A  and  B  at  t-5.  Beyond  a  radial  distance  of  0.5  inches,  the 
material  is  void- free.  The  strength  relaxation  can  be  seen  between  A'  and  B' 
due  to  this  small  amount  of  void  generation.  As  the  void  distribution  spreads 
to  larger  distances  (between  C  and  D)  by  20  microseconds,  the  material 
strength  diminishes,  as  represented  by  the  solid  line  between  C'  and  D'  over 
larger  radius . 
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Figure  40.  The  Loading  Path  (solid  line)  with  Respect  to  the  Pres 
Dependent  Yield  Surfaces. 
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The  void  volume  fraction  contours  for  different  time  intervals  are 
plotted  in  Figure  42.  At  5  microseconds,  void  nucleation  and  growth  occurred 
at  a  distance  0.9  inches  from  the  top  of  the  target.  Since  the  material 
around  this  region  experiences  relatively  large  tensile  stresses,  the  spall 
initiates  in  this  region  first.  The  growth  rate  was  relatively  lower  outside 
the  region  enclosed  by  the  contour  in  Figure  42a.  Between  5  and  10 
microseconds,  the  void  growth  extended  to  a  larger  region  as  shown  in  Figure 
42b.  At  15  microseconds,  the  spall  process  is  completed  and  intense  void  con¬ 
centration  occurred.  The  inner  contours  in  Figure  42c  represent  this  intense 
spall  zone.  The  void  volume  fraction  in  this  zone  was  well  above  0.5.  In  a 
penetration  experiment,  the  material  often  physically  separates  and  forms 
stress  free  fracture  surfaces  in  the  target  due  to  void  coalescence;  in  some 
cases,  spall  fragments  are  ejected  out  from  the  back  of  the  target. 

For  completion,  the  case  of  a  penetration  process  without  any  material 
degradation  was  simulated.  In  the  calculation,  the  void  volume  fraction  was 
set  to  zero  and  also  no  other  spall  criterion  was  included.  The  material  be¬ 
havior  was  therefore  described  by  the  JC  or  BP  viscoplastic  models  and  without 
the  RDG  model.  In  Figure  43,  the  matrix  and  aggregate  strength  history  plots 
from  the  two  analyses,  one  with  material  degradation  (RDG  failure  model)  and 
the  other  without  any  degradation  model,  are  presented.  Again  these  plots  are 
for  the  element  #736  of  the  finite  element  mesh  shown  earlier  in  Figure  37. 
Since  the  aggregate  and  matrix  strengths  are  the  same  in  the  absence  of  any 
porosity,  the  strength  remained  high  as  represented  by  the  dashed  line  in 
Figure  43. 

The  deformed  configurations  of  the  penetrator  and  the  target  at  time  20 
microseconds  are  compared  in  Figure  44.  The  target  exhibited  a  relatively 
larger  bulge  from  enhanced  plastic  flow  due  to  void  softening  of  the  material 
Since  the  target  material  in  front  of  the  projectile  has  been  damaged  due  to 
spalling,  the  penetration  process  is  accelerated  and  this  allows  the  projec¬ 
tile  to  penetrate  easier.  When  there  is  no  material  degradation  due  to  void 
nucleation  and  growth,  the  target  material  is  relatively  stronger,  and  there¬ 
fore  the  penetration  rate  is  slower  compared  to  the  case  in  which  the  material 
degrades.  The  differences  in  the  deformed  target  shapes  and  the  penetration 
rates  demonstrate  the  importance  of  accurate  dynamic  failure  model. 
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One  could  expect  dependable  and  realistic  failure  predictions  from 
physically  based  evolution  equations  of  the  failure  model.  Since  the  damage 
evolution  depends  on  the  time  and  loading  history,  use  of  appropriate  damage 
evolution  laws  and  viscoplastic  description  of  the  plastic  flow  are  extremely 
important . 

6.3  Spall  in  a  Solid  Cone  Target 

Bless  [41]  considered  a  plate  impact  configuration  in  which  a  flyer 
plate  impacts  the  base  of  a  solid  right  circular  cone  as  schematically  shown 
in  Figure  45a.  The  target  was  a  cone  of  25  mm  in  height  with  a  90°  included 
angle.  The  flyer  and  the  target  materials  were  1020  steel.  The  flyer  thick¬ 
ness  and  diameter  were  2  mm  and  38  mm  respectively.  Depending  on  the  impact 
velocity,  a  complex  spall  pattern  was  generated  in  the  cone.  The  pattern  con 
sisted  of  tunnel  cracks  parallel  to  the  axis  and  longitudinal  cracks  parallel 
to  the  circumference  of  the  cone,  as  shown  in  Figure  45b.  This  cone  impact 
configuration  is  a  two-dimensional  axisymmetric  geometry.  However,  the  crea¬ 
tion  of  spall  regions  leading  to  macrocrack  formations  destroy  the 
ax isymme tricity ,  and  the  configuration  eventually  becomes  three  dimensional. 
The  EPIC-2  code  analysis  includes  only  the  onset  of  spallation  and  maps  the 
various  spall  zones  in  the  target.  A  flyer  piate  of  diameter  50  mm  was  as¬ 
sumed  in  the  simulation.  The  high  strain  rate  behavior  of  1020  steel  was 
modeled  using  the  Bodner- Partom  model.  The  corresponding  model  constants  are 
tabulated  earlier  in  Table  1.  The  RDG  model  constants  have  been  determined 
from  plate  impact  experimental  data  and  are  given  in  Table  4.  The  solution 
was  carried  out  for  60  microseconds.  The  finite  element  mesh  for  the  cone 
configuration  at  t-0  and  t-60  microseconds  are  shown  in  Figure  46. 

EPIC  analysis  showed  formation  of  incipient  spall  surfaces  within  five 
microseconds  after  the  impact  due  to  voids  nucleation  and  growth.  Figure  47 
shows  the  void  volume  fraction  contours  for  t-5  and  t-60  microseconds. 
Comparisons  between  these  two  time  intervals  clearly  indicated  that  the  spall 
evolution  is  almost  complete  during  the  initial  five  microseconds.  A  qualita 
tive  comparison  between  the  simulation  and  the  experiment  (Figure  45b)  clearl 
shows  that  the  RDG  model  could  reproduce  the  salient  features  of  the  spall 
fracture  evolution. 
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Figure  45.  (a)  Schematic  Two  Dimensional  Sketch  of  a  Flyer  Plate  Impact  on 

the  Base  of  a  Right  Solid  Cone.  (b)  Recovered  and  Sectioned  Solid 
Cone  Revealing  the  Complex  Spall  Zones. 


88 


for  the  Cone  Impact;  Geometry  at  ( a ) 


7.  SUMMARY  AND  CONCLUSIONS 


The  RDG  model  is  a  continuum  mechanics  based  ductile  failure  model. 

The  model  formulation  is  three  dimensional  and  therefore  applicable  under  a 
general  stress-strain  state.  The  RDG  model  describes  several  fundamental 
aspects  of  the  failure  process.  When  the  failure  process  is  initiated  due  to 
void  nucleation  in  the  material,  the  plastic  flow  becomes  pressure  dependent. 
Also,  under  high  velocity  impact  loading,  most  metals  exhibit  strain  rate  de¬ 
pendent  behavior.  The  RDG  model  includes  the  effects  of  pressure  and  strain 
rate  in  the  constitutive  (strength)  model  formulation.  An  important  feature 
of  the  RDG  model  is  its  accurate  modeling  of  the  intact  material  using  ap¬ 
propriate  viscoplastic  equations.  The  Bodner- Partom  (BP)  equations  were 
selected  for  this  purpose.  The  state  variables  in  the  RDG  model  formulation 
become  the  equivalent  plastic  strain  rate  in  the  matrix  material  and  the  BP 
model  state  variable  Z,  which  describes  the  loading  history  effects. 

A  pressure  dependent  yield  function  serves  as  the  plastic  potential  in 
the  derivation  of  plastic  strain  rate  for  the  void  contained  aggregate.  Since 
the  aggregate  plastic  strain  rates  are  directly  related  to  the  equivalent 
plastic  strain  rate  of  the  matrix  material,  the  strain  rate  and  loading  his¬ 
tory  effects  enter  into  the  failure  model  formulation  through -the  Bodner - 
Partom  model.  The  RDG  failure  model  does  not  require  separate  loading  or 
unloading  conditions.  Both  elastic  and  plastic  rate  components  are  taken  to 
be  always  nonzero,  and  the  same  relations  are  intended  to  hold  under  loading 
and  unloading  conditions  as  in  the  Bodner- Partom  model.  When  the  equivalent 
plastic  strain  rates  are  numerically  very  small  ( <0.0001),  the  deformation  is 
considered  to  be  elastic. 

Another  appealing  feature  of  the  RDG  failure  model  is  the  void  growth 
law.  Assuming  that  the  void  volume  fraction  rate  is  equal  to  the  summation  of 
the  normal  plastic  strain  rate  components  (dilation)  removes  the  requirement 
for  a  new  evolution  law  for  the  void  growth  rate.  This  modeling  approach  sig¬ 
nificantly  reduces  the  number  of  model  constants.  While  the  Bodner-Partom 
model  requires  a  number  of  constants,  the  determination  of  those  constants  is 
fairly  straight  forward  and  Rajendran  et  al.,  [24]  reported  a  procedure  to 
calibrate  the  constants  using  data  from  a  limited  number  of  standard  impact 
tests . 
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In  the  RDG  model,  there  are  three  constants  for  the  nucleation,  and  two 
constants  for  the  growth  of  voids.  Model  constants  were  determined  from  the 
numerical  simulations  of  a  plate  impact  experiment.  The  standard  deviation  is 
arbitrarily  set  to  one  third  or  one  fourth  of  this  stress.  The  initial  value 
for  the  void  volume  fraction  and  the  nucleation  threshold  stress  can  be  ad¬ 
justed  to  match  the  time  arrival  of  the  spall  signal.  The  initial  slope  and 
amplitude  of  the  spall  signal  will  guide  the  selection  of  void  growth  con¬ 
stants  . 

The  RDG  model  was  implemented  into  the  EPIC-2  finite  element  code.  A 
numerically  stable  and  accurate  solution  technique  based  on  a  diagonally  im¬ 
plicit  Runge-Kutta  (DIRK)  method  was  used  to  solve  the  complex  and  stiff 
governing  differentia]  equations.  Various  aspects  of  the  numerical  scheme 
were  demonstrated  in  References  1  and  20.  For  solving  only  the  Bodner- Partom 
equations,  an  iterative  radial  return  scheme  has  been  successfully  imple¬ 
mented.  The  BP  model  constants  had  already  been  determined  for  these 
materials  and  reported  by  Rajendran  and  Bless  [41],  The  model  accurately 
matched  several  of  the  plate  impact  experiments.  The  RDG  model  constants  for 
OFHC  copper,  tantalum,  Armco  iron,  and  1020,  MAR  200,  MAR  250,  AF1410 ,  C1008, 
and  HY100  steels  were  determined  and  presented  in  the  report. 

The  capability  of  the  RDG  model  was  extended  to  multiple  impact  loading 
conditions.  The  model  could  accurately  predict  the  failure  process  in  a 
twice  -  impacted  target  plate.  Both  void  growth  and  collapse  were  modeled  by 
the  RDG  model.  To  describe  the  pore  collapse  correctly,  the  critical  value  of 
the  void  volume  fraction  was  calibrated  for  collapse.  For  complete  spalla¬ 
tion,  the  critical  void  volume  fraction  was  arbitrarily  assumed  to  be  1.0  for 
most  metals. 

The  results  from  the  analysis  of  application  problems  require  addi¬ 
tional  in-depth  analyses  supported  by  experimental  data.  For  example,  to 
substantiate  the  tensile  necking  analysis,  high  speed  photographic  split 
Hopkinson  bar  tests  on  a  shallow  notch  specimens  should  be  conducted.  The  ex¬ 
perimentally  determined  neck  profiles  should  be  compared  with  the  RDG  model 
predictions.  For  this  purpose,  selection  of  a  material  in  which  microvoids 
have  been  observed  under  dynamic  tensile  loading  conditions  is  recommended. 
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It  is  desirable  to  perform  instrumented  penetration  experiments  for 
validating  the  RDG  model  prediction.  Radiographs  of  the  projectile  and  target 
during  the  penetration,  and  post- failure  analysis  on  the  target  would  identify 
additional  model  requirements.  These  radiographs  are  particularly  important 
for  verifying  the  spall  fragmentation  studies. 
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APPENDIX  A 

Numerical  Solution  of  the  Governing  Equations  in  Section  3 

The  governing  equations  described  in  Section  3  are  rearranged  and 
numerically  integrated.  For  stability,  accuracy,  and  uniqueness  of  the  solu¬ 
tion,  a  second  order  diagonally  implicit  Runge-Kutta  (DIRK)  method  was 
considered  for  integrating  the  system  of  stiff  ordinary  differential  equa¬ 
tions.  The  DIRK  scheme  is  outlined  in  Appendix  B.  The  DIRK  variables  are 
integrated  using  the  estimated  rates  for  the  current  time.  Grove  et  al.  [38] 
successfully  implemented  the  RDG  model  in  the  EPIC-2  code  using  the  DIRK 
scheme . 

The  EPIC-2  code  calculates  the  current  strain  rates  in  each  material 
element.  Before  any  void  nucleatfon  occurs,  the  volumetric  strain  of  the 
matrix  and  the  aggregates  are  the  same.  However,  once  void  nucleation  occurs 

in  an  element,  the  incremental  volumetric  strain  of  the  aggregate,  e  is  the 

0 

sum  of  the  incremental  elastic  volume  change,  ,  in  the  aggregate  (in  the 
same  sense  as  Mackenzie's)  and  the  incremental  volume  change  due  to  the  vis¬ 
coplastic  growth  of  voids,  .  This  relationship  can  be  related  directly  to 

tne  aggregate  elastic  compressible  strain,  n  ,  from  Equation  (27)  as  follows. 

a6 

The  relationship  between  the  true  volumetric  strain  and  the  compressible 
strain  is  by  definition, 

1  +  n  -  (1  -  fQ)  exp(  -«v)  .  (A. 1) 

where  f^is  the  initial  void  volume  fraction  at  the  beginning  of  the  simula¬ 
tion.  The  relationship  between  the  true  plastic  dilatation  and  the  void 
volume  fraction  is  given  by  the  direct  integration  of  Equation  (10)  as, 


1  -  f  -  (1  -  fQ)  exp(  - e?^  ) 


(A. 2) 


where  the  plastic  dilatation,  ,  is  zero  at  f  -  fQ  .  Substituting  equations 
( A . 1 )  and  (A. 2)  into  Equation  (27),  we  obtain  the  relationships, 


u  -  exp(  -€  )  -  1 
ag  v 

e  p 

e  -  e  -  £., 

V  V  u 


and 


(A. 3) 
(A. 4) 
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£®  is  the  true  elastic  dilatation  strain  of  the  aggregate.  This  development 

*  p 

suggests  that  the  plastic  strain  rates,  «^p  as  given  by  Equation  (20)  should 
be  numerically  integrated  instead  of  integrating  Equation  (23).  The  aggregate 
pressure,  P,  is  evaluated  by  using  Equation  (A. 3)  and  the  procedure  outlined 
in  Sections  3.4  and  3.5. 


The  true  elastic  deviatoric  strains  of  the  aggregate  are  given  by  the 
equation , 


(« 


ki/3>hj  1 


(A. 5) 


This  leads  to  the  degraded  deviatoric  stresses  as  in  the  equation, 


S.  . 


2  G  e 


ij 


(A. 6) 


One  numerical  bonus  of  this  approach  is  that  the  use  of  Equation  (A. 2)  to  cal¬ 
culate  f  avoids  numerical  integration  of  the  stiff  differential  equation  (10). 
Equations  (2)  and  (7)  lead  to  stiff  differential  equations;  therefore,  the 
DIRK  scheme  permits  extremely  small  time  steps  to  maintain  solution  stability 
and  accuracy.  For  numerical  efficiency  (i.e.  to  avoid  extremely  small  time 
steps),  Equations  (2)  and  (7)  were  analytically  integrated  to  result  in  the 
following  equations  for  the  BP  model  loading  history  parameter  and  the  void 
nucleation  volume  fraction, 


Z  -  -  (Z^-  Zq)  exp (  -terra  )  (A. 7) 

term  -  m^W^  +  (ra^/o)  [l  -  exp(  -a  W  )  ]  (A  8) 

fn  -  (fj/2)  [  erf  ( (P  +Ym-  s^)  -  erf  ((-oN)/(/2  S]_))]  + 

(f2/2)  [ erf ( (D^-  eN)/(7?  s2))  -  erf  ((-eN)/(7?  s2))]  (A. 9) 

One  numerical  bonus  of  having  Equation  (A. 9)  is  not  needing  to  evaluate  the 
effective  matrix  stress  rate,  Y  ,  as  was  required  in  Equation  (7).  The  value 
for  the  effective  matrix  stress  is  obtained  from  the  yield  function,  Equation 
(16),  as  in  the  equation, 

2 

Y*  -  [  (2+p2)  J2  +  I*]/  5(p)  (A. 10) 
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The  yield  function  described  by  Equation  (16)  was  treated  as  a  plastic  poten¬ 
tial.  Implicitly  assumed  was  that  the  loading  and  unloading  are  controlled  by 
the  BP  model  through  the  matrix  equivalent  plastic  strain. 


This  rearrangement  of  the  governing  equations  has  resulted  in  a  fewer 

differential  equations  which  are  relatively  less  stiff.  Thus,  the  numerical 

integration  problem  was  reduced  to  the  solution  of  the  rate  equations  for  , 

ra 

W  ,  I ,  and  «?,  . 

P  ij 

•  p 

Examination  of  Equation  (1)  shows  that  is  ultimately  the  function  of 

m 

the  same  four  integration  variables  (the  variables  input  from  EPIC  are  con¬ 
sidered  constant  during  the  EPIC  time  step).  That  is,  the  Z  in  Equation  (1) 

is  onlv  a  function  of  W  through  Eqs .  (A. 7)  and  (A. 8).  The  Y  in  Equation  (1) 

p  m 

is  a  function  of  f,  P,  and  S  through  Equation  (A. 10).  The  variables  f,  P, 

and  S„  are  in  turn  functions  of  through  Eqs.  (24)  to  (26)  and  Eqs  (A. 2) 

to  (A. 6) .  Since  the  variables  P  and  S, .are  also  direct  functions  of  f  in  the 

ij  n 

degraded  modulus  and  that  f^  is  in  turn  a  function  of  P  and  through 
Equation  (A. 9),  an  efficient  convergent  iterative  scheme  was  developed  to 
solve  f^  simultaneously  with  P  and  .  The  temperature  in  Equation  (4)  and 
the  pressure  in  Equation  (26)  are  also  direct  functions  of  the  specific  inter¬ 
nal  energy.  The  rate  equation  for  the  internal  energy  is, 


I  -  (-P  £v  +  S.^.jVUnO  (A. 11) 

Examination  of  these  equations,  Equation  C18)  for  the  plastic  work 
rate  and  Equation  (20)  for  the  aggregate  plastic  strain  rates,  show  that  the.’ 
are  also  ultimately  functions  of  the  same  four  dependent  integration  vari¬ 
ables,  D*5  ,  W  ,  I  ,  and  £?.. 

m  p  lj 

The  corresponding  rate  equations  were  solved  using  the  DIRK  scheme. 
Demonstrated  in  Section  3,  was  the  solution  stability,  accuracy,  and  unique¬ 
ness  with  respect  to  different  finite  element  mesh  size  and  time  step  through 
one  dimensional  plate  impact  solutions.  This  was  again  demonstrated  through 
two  dimensional  solutions  of  several  application  problems. 
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APPENDIX  B 

The  Diagonally  Implicit  Runge  Kutta  (DIRK) 
Scheme  With  Time  Step  Control 


31 ; 


The  diagonally  implicit  Runge-Kutta  integration  is  given  by  the  formula 


K(.n)-  f(t  +  c.h  ,  y  +  h  S  a..K[n)).  i-1  ( 1 ) q 

i  n  l  7n  J  lj  j  M 


Yn+1 


-  y  +  h  S  b .  K. 
J  n  Til 


(n) 


(dv/dt)  -  f(t  ,  y  ) 
n  n  n 


(B.  1) 
(B.2) 
(B.3) 


h  is  the  time  step  and  y^  is  the  vector  of  the  integrated  variables  at  time, 

t  Note  that  K.  is  also  a  function  of  itself,  which  requires  an  iterative 
n  i  n 

corrector  procedure.  That  is,  the  calculated  value  of  K.  is  re  -  subs t i tuted 

l 

into  the  nonlinear  equation  (B.l)  until  convergence  is  achieved  within  an  er¬ 
ror  tolerance.  If  convergence  is  net  achieved  after  a  certain  number  of  re- 
substitutions  of  ,  then  the  time  step,  h,  is  reduced.  Convergence  of  K.  is 
assured  if  a  small  enough  time  step  is  chosen.  In  the  RDG  model  implementa¬ 
tion,  if  convergence  was  achieved  within  a  very  small  error  tolerance,  the 
time  step  h  was  gradually  increased  to  improve  computation  time.  This  scheme 
was  implemented  into  the  shock  wave  propagation,  finite  element  code,  EPIC- 2 


The  EPIC-2  host  code  is  based  on  the  time  -  centered  integration  approach 
with  the  time  step  controlled  by  the  Courant  stability  criterion;  therefore, 
the  time  steps  usually  remain  small.  Thus,  the  lowest  order  DIRK,  (q-1),  will 
be  adequate  for  stable  and  accurate  solutions.  For  the  values,  q-1,  a^-1/2  , 
b^-1,  and  c^-1/2,  the  DIRK  method  becomes  second-order  accurate.  This  special 
case  is  equivalent  to  Bass  and  Oden's  [B2]  recommendation  of  their  Euler 
predictor  -  trapezoidal  corrector  scheme.  Saliba  [B3]  has  effectively  proved 
this  special  ca.e  as  stable  and  accurate  in  conjunction  with  the  mechanical 
equilibrium  equations  for  quasi-static  problems. 

With  the  DIRK  method  applied  to  the  ODE  of  any  viscoplastic  formula¬ 
tion,  a  single  corrector  step  was  found  sufficient  for  accurate  calculations 
most  of  the  time.  However,  during  the  stiff  phase  of  the  solutions,  such  as 
the  unloading,  the  DIRK  time  steps  require  significant  reduction  even  with 
several  re  -  subs t i tut  ions  of  K  Maximum  efficiency  during  the  stiff  phases 
was  achieved  as  follows.  If  after  three  re-substitutions  of  K^  the  relative 
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error  is  not  within  two  percent,  the  DIRK  time  step  is  reduced  by  one-tenth. 

The  process  continues  until  the  two  percent  error  is  obtained.  If,  however, 

after  three  re  -  substitutions  of  K.  the  relative  error  is  within  two-tenth  of 

i 

a  percent,  the  DIRK  time  step  is  increased  in  an  arbitrary  manner  by  332  for 
the  next  time  step.  This  process  continues  until  the  end  of  the  EPIC  time 
step  is  reached. 
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